DISC  FSLE  COPY 


Computing 
\  About 
\S\  Physical 


DT1C 

8ELECTE 

NQV28199C 


I 


— y — y  -V- — y — V 

\  \  \  \  \  \  \  \  x 


\  \ \ V  \ \ \ \  \ \ \  \  \ \ \ \\\\  \ \ \  \ \ \ 

\  \  .  ■ ,  \  \  \  \  ,  \  \  \  \  \  \  \  \  \  \  \  \  \  \  \  \ 

■>-  V-\  -V \  \-\  \-v-\-V-V  v\\\\  -\  VV\  \\  \  - 

x  \  \  x.  \  \  \  X.  \  \  \  \  \\  \\  \  \  \  \  \  \  \  \  X 


\  \  •  'x  \  \  \  \  \  \ 
\  \  \ \ \\ \  \\  • ' 


Purdue  ( ’nil  vrsity 
Department  of 
\  \  Computer  Sciences 


,  \  \  x  ""  \  \  \\  \\  \\\  \\\  \  \  \\  \  \\\ 

v  ^  v\m\ w\\^\\\x\v\ 

x  -  X-  \  V  \  X  '  '•  V  X  \  A-  \  \  \  - \  \  \  \  -  \ -  \-v  x  \  X  - 

■  \  \  \  \  \  ^  \\ \  \\W  \  v  \  \ \  \\  \  \  x 


•Apsxa—A  for  pnbfic  wbail 
K4dbaet»  Utt&Tai!*d 


8 


SURFACE  APPROXIMATIONS 
IN  GEOMETRIC  MODELING* 


Jung-Hong  Chuang 

Computer  Sciences  Department 
Purdue  University 
Technical  Report  CSD-TR-1023 
CAPO  Report  CER-90-37 
September,  1990 


1 


*  A  Thesis  submitted  to  the  faculty  of  Purdue  University  in  partial  fulfillment  of  the  require¬ 
ments  for  the  degree  of  Doctor  of  Philosophy,  August  1990. 


SURFACE  APPROXIMATIONS  IN  GEOMETRIC  MODELING 


A  Thesis 

Submitted  to  the  Faculty 
of 

Purdue  University 
by 


Jung-Hong  Chuang 


In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree 

of 

Doctor  of  Philosophy 


A  i  must  1900 


11 


To  my  parents 


ACKNOWLEDGMENTS 


I  would  like  to  express  my  deepest  gratitude  to  my  advisor,  Professor  Christoph 
M.  Hoffmann,  for  his  guidance  and  encouragement,  and  for  his  detailed  reading 
and  valuable  feedback.  He  has  significantly  influenced  my  professional  conduct 
and  training.  Professors  Robert  Lynch.  Bradley  J.  Lucier  and  Elias  N.  Houstis 
deserve  special  thanks  for  their  interest  in  my  research  and  for  serving  on  my 
thesis  committee.  I  am  also  thankful  to  Professor  Chanderjit  Bajaj  for  his  interest 
with  several  aspects  of  my  study.  Thanks  to  Kevin  Kuehl  for  his  great  effort 
to  implement  the  algorithms  on  piecewise  linear  approximations.  Thanks  also  to 


Pamela  Vermeer  for  reviewing  this  manuscript. 

I  thank  all  the  friends  who  have  offered  assistance  and  encouragement  during 
my  stay  at  Purdue.  Notable  among  these  are  Danny  Chen,  Ching-Shoei  Chiang, 
Shy-Renn  Lian,  Hsin  Pan.  Vic  Pollara.  Preeti  Shrikhande,  Charles  and  Jui-Fen 
Trappey,  Jyh-Jong  Tsay,  Ko-Yang  Wang,  JiaXun  Yu,  and  Jianhua  Zhou. 

My  wife  Whey- Ming  Tina  Song  deserves  special  thanks  for  her  encouragement, 
love,  patience,  and  understanding.  Such  as  they  are.  my  accomplishments  belong 
equally  to  her. 

Finally,  to  my  parents,  who  seemed  to  be  proud  of  what  I  did.  I  would  like  to 
express  my  deepest  appreciation.  This  thesis  is  dedicated  to  them. 

I  would  like  to  acknowledge  the  financial  support  provided  by  the  Office  of 


Naval  Research  under  contracts  N00014-S6-K-0465  and  N000014-90-J-1599.  and 


the  National  Science  Foundation  under  grant  CCR-S6-19S17. 


Statement 

Tilborg. 

1133. 


VHG 


11/28/90 


iv 


TABLE  OF  CONTENTS 

Page 

LIST  OF  FIGURES .  vi 

LIST  OF  TABLES . viii 

ABSTRACT .  ix 

1.  INTRODUCTION .  1 

1.1  Surface  Representations  in  Solid  Modeling .  2 

1.1.1  Parametric  Surfaces .  3 

1.1.2  Implicit  Surfaces .  3 

1.1.3  Implicit  Surfaces  in  High-Dimensional  Space .  4 

1.2  Implicit  Approximations  of  Parametric  Curves  and  Surfaces  ....  10 

1.3  Local  Approximations  of  2-D  Surfaces .  13 

1.4  Piecewise  Linear  Approximations  of  2-D  Surfaces .  13 

1.5  Some  Preliminaries .  IS 

1.6  Thesis  Organization .  20 

2.  IMPLICIT  APPROXIMATION  OF  CURVES  AND  SURFACES  21 

2.1  Local  Implicit  Approximation  of  Parametric  Plane  Curves .  21 

2.1.1  A  Recurrence  for  a*,.  23 

2.1.2  Derivation  of  the  Method  .  25 

2.1.3  Error  Analysis .  32 

2.1.4  Experiments .  35 

2.2  Local  Implicit  Approximation  of  Parametric  Surfaces .  44 

2.2.1  A  Recurrence  for  a,; .  4S 

2.2.2  Derivation  of  the  Method  .  50 

2.3  Remarks  on  Resultants .  53 


G«  G<  Gi 


V 


Page 

3.  LOCAL  APPROXIMATIONS  OF  2-D  SURFACES  .  57 

3.1  Local  Geometry  of  the  Projected  Surface .  58 

3.1.1  Normal  Vector  .  59 

3.1.2  Tangent  Vectors .  60 

3.1.3  Normal  Curvature  .  64 

3.2  Degree  Two  Implicit  Approximation .  70 

3.3  Local  Explicit  Approximation .  74 

3.4  Local  Parametric  Approximation .  77 

4.  PIECEWISE  APPROXIMATIONS  OF  2-D  SURFACES  ....  82 

4.1  The  PL  A  of  Offset,  Voronoi  and  Blending  Surfaces  .  S3 

4.2  Proposed  Approach .  86 

4.3  Detailed  Computations .  89 

4.3.1  Intersecting  A  Parametric  Approximant  with  Edges  of  a  Cube  S9 

4.3.2  Intersecting  A  Projected  Surface  with  Edges  of  a  Cube  ...  94 

4.3.3  Strategies  for  Stepping .  99 

4.3.4  Newton  Iteration . 102 

4.3.5  Some  Error  Analyses . 103 

4.3.6  Adaptive  Subdivision . 104 

4.3.7  Polygonization  and  Local  Refinement . 107 

5.  CONCLUSIONS  AND  FUTURE  WORK . 109 

.1  Local  Implicit  Approximations  of  Curves  and  Surfaces . 109 

.2  Local  Approximations  of  2-Surfaces . 110 

.3  Piecewise  Approximations  of  2-Surfaces . 110 

BIBLIOGRAPHY . 112 


VITA 


118 


VI 


LIST  OF  FIGURES 

Figure  Page 

2.1  Error  estimate .  34 

2.2  Curve  c4(t)  .  37 

2.3  Curve  c5(t)  .  38 

2.4  Curve  e6(t)  .  39 

2.5  Curve  c7(t)  .  40 

2.6  Curve  cs(t)  .  45 

2.7  Curve  c9(t)  .  46 

2.8  Curve  c9(t)  .  47 

2.9  /4  =  0  n  h  =  0  and  g2  =  0  D  h  =  0 .  54 

2.10  /4  =  0  n  h  =  0  and  g3  =  0  n  h  =  0 .  55 

4.1  Intersecting  <&i(s.t)  with  a  cube  B  (Regular  cases) .  90 

4.2  Intersecting  $i(s,<)  with  a  cube  B  (Degenerate  cases) .  91 

4.3  Intersecting  a  tangent  plane  with  the  cube  B .  92 

4.4  An  approximate  point  is  refined  to  a  surface  point .  95 

4.5  An  approximate  point  is  refined  to  two  surface  points .  96 

4.6  Two  approximate  points  are  refined  to  one  surface  point .  96 

4.7  The  refinement  of  t, i ) .  98 

4.8  Special  cases  for  stepping . 101 


vii 

Figure  Page 

4.9  An  invalid  polygon . 105 

4.10  A  crack  on  the  shared  face . 106 


via 


LIST  OF  TABLES 


Table  Page 

2.1  Error  of  approximation  for  curve  c^{t)  41 

2.2  Error  of  approximation  for  curve  c $(t)  42 

2.3  Error  of  approximation  for  curve  C6(£)  43 

2.4  Maximal  deviations  between  the  intersection  curves .  53 


ABSTRACT 


Chuang,  Jung-Hong.  Ph.D.,  Purdue  University,  August  1990.  Surface  Approxi¬ 
mations  in  Geometric  Modeling.  Major  Professor:  Christoph  M.  Hoffmann. 

One  of  the  major  research  efforts  in  the  field  of  solid  modeling  focuses  on  ex¬ 
tending  the  geometric  coverage  of  modeling  systems  and  on  incorporating  complex 
free-form  surfaces.  Some  major  obstacles  to  this  goal  include  computing  and  rep¬ 
resenting  intersection  curves  of  two  general  surfaces,  and  computing  and  rendering 
very  complex  surfaces,  including  offset,  Voronoi.  and  blending  surfaces.  We  present 
local  and  global  approximation  schemes  that  are  expected  to  be  of  practical  value 
in  overcoming  the  above  problems.  For  parametric  curves  and  surfaces,  we  present 
a  method  for  computing  an  implicit  approximant  of  low  degree  that  approximates 
the  curves  or  surface  locally  and  achieves  an  order  of  contact  that  can  be  prescribed 
in  advance.  In  principle,  the  method  is  capable  of  exact  implicitization.  ^Several  \ 
surfaces,  including  offsets,  blends,  and  Voronoi  surfaces  can  be  defined  as  the  nat¬ 
ural  projections  to  R3  of  2-surfaces  in  R”.  n  >  3.  The  2-surface  in  R’4  is  the  zero 
set  of  a  system  of  nonlinear  equations  in  n  variables.  We  present  algorithms  that 
compute  the  normal,  tangent  vectors,  and  normal  curvatures  of  the  projected  sur¬ 
face  directly  from  the  nonlinear  system  without  variable  elimination.  Methods  are 
presented  as  well  that  compute  the  explicit  and  parametric  approximations  of  the 
projected  surface  locally.  Finally,  for  a  given  2-surface  in  Rn,  n  >  3,  an  algorithm 
is  given  that  computes  the  piecewise  linear  approximation  of  the  projected  surface 
globally  with  all  major  computations  performed  in  3-space. 
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1.  INTRODUCTION 

Solid  modeling  has  received  much  attention  throughout  the  academic  and  in¬ 
dustrial  communities  for  nearly  three  decades.  Even  though  significant  progress 
has  been  made  in  basic  research  and  in  the  capabilities  of  commercially  available 
solid  modelers,  many  current  solid  modeling  systems  allow  only  severely  geomet¬ 
ric  primitives.  For  example  PADL,  a  solid  modeler  developed  at  the  University 
of  Rochester  by  Requicha  and  Voelcker.  only  allows  planar,  spherical,  cylindrical, 
conical,  and  toroidal  faces  [55]. 

Recent  solid  modeling  research  efforts  are  being  directed  at  extending  the  ge¬ 
ometric  coverage  of  solid  modelers  and  incorporating  complex  free-form  surfaces. 
Some  major  obstacles  to  this  goal  have  been  computing  and  representing  intersec¬ 
tion  curves  of  two  general  surfaces,  and  computing  and  rendering  very  complex 
surfaces,  including  offset,  Voronoi,  and  blending  surfaces.  In  this  thesis,  local  and 
global  approximation  schemes  are  presented  that  are  expected  to  be  of  practical 
value  in  overcoming  the  problems  of  geometric  coverage. 

An  algebraic  surface  in  three  dimensional  space  can  be  given  by  an  implicit 
representation  as  the  polynomial  equation  g{x,  y,z)  =  0.  Some  algebraic  surfaces 
can  also  be  given  parametrically  as  x  =  hi{s,t),  y  =  h2{s,t),  :  =  A3(s,<),  where 
the  /i,  are  polynomials  or  ratios  of  polynomials.  Since  both  representations  have 
their  own  strengths  and  weaknesses,  many  geometric  computations  could  become 
simpler  or  practical  if  both  representations  were  available  and  their  complemen¬ 
tary  strengths  could  be  utilized.  Thus,  the  problem  of  how  to  convert  from  one 
representation  to  the  other  is  of  great  practical  importance.  However,  conversions 
between  them  are  not  always  possible.  While  general  techniques  exist  for  con¬ 
verting  from  parametric  to  implicit  form,  by  a  process  called  imphcitization  based 
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on  classical  elimination  theory,  e.g.,  [46,  47,  49,  5S,  65],  only  a  subset  of  implicit 
curves  and  surfaces  have  a  parametric  representation.  Moreover,  the  techniques 
for  implicitization  are  extremely  expensive  and  hence  their  practical  use  is  quite 
limited.  Since  exact  conversions  between  representations  might  be  impossible  or 
impractical,  local  approximations  in  terms  of  parametric  or  implicit  forms  should 
be  valuable.  We  have  derived  several  local  approximation  techniques  for  curves 
and  surfaces  that  are  defined  parametrically  or  implicitly,  and  expect  that  these 
techniques  will  be  of  practical  interest  when  they  are  incorporated  into  algorithms 
for  surface  interrogation  such  as  computing  surface  intersections  [23,  24], 

Offset.  Voronoi,  and  blending  surfaces  can  be  constructed,  with  an  algebraic  for¬ 
mulation.  as  the  natural  projection  to  R3  of  2-dimensional  manifolds  (2-surfaces)  in 
a  higher-dimensional  space.  For  visualization  purposes,  a  piecewise  linear  approx¬ 
imation  (PLA)  of  these  surfaces  is  desirable.  As  a  global  approximation,  the  PLA 
of  parametric  surfaces  has  been  extensively  studied  and  is  a  widely  used  tool  for 
rendering  surfaces,  as  well  as  for  evaluating  surface  intersections.  Comparatively, 
much  less  attention  has  been  paid  in  the  literature  to  piecewise  approximations  of 
implicitly  defined  surfaces  [7,  17,  57].  We  present  new  algorithms  for  computing 
PLAs  of  surfaces  defined  implicitly  by  sets  of  equations,  including  offset  surfaces, 
Voronoi  surfaces,  and  spherical  blending  surfaces. 

1.1  Surface  Representations  in  Solid  Modeling 

2-surfaces  can  be  represented  either  in  parametric  or  in  implicit  form.  The 
implicit  representation  is  a'  .active  for  determining  directly  whether  a  point  is  on 
the  surface  by  checking  whether  or  not  it  sat.sfies  the  implicit  equation.  For  the 
parametric  representation,  on  the  other  hand,  it  is  much  easier  to  generate  points 
on  the  surface.  Major  parametric  surfaces.  includinQ  Bezier  surfaces,  and  nonuni- 
forrn  rational  B-splines.  are  attractive  in  interactive  design  because  the  manner 


in  which  coefficient  changes  alter  the  surface  shape  qualitatively  can  be  grasped 
intuitively.  However,  this  is  presently  not  the  case  for  implicit  surfaces. 

1.1.1  Parametric  Surfaces 

The  parametric  representation  of  a  2-surface  in  R3  is 

x  =  x(u,u) 

U  =  y(u,  v)  (1.1) 

-  =  ;(n.y) 

and  is  the  range  of  a  map  from  R2  to  R3,  where  u  and  v  are  usually  restricted 
to  a  standard  domain,  say  to  the  unit  square  [0,1]  x  [0.1].  These  parameter 
limits  define  a  bounded  rectangular  piece  of  surface,  or  a  surface  patch.  The 
functions  x(u.v),  y(u,v),  and  z(u,v)  customarily  are  polynomials  or  ratios  of 
polynomials  in  u  and  v.  For  the  major  parametric  surfaces  used  in  Computer 
Aided  Geometric  Design  (CAGD),  the  polynomials  are  represented  in  a  particular 
basis,  for  instance,  in  the  Bernstein-Bezier  basis.  This  allows  us  to  relate  the 
coefficients  of  the  coordinate  functions  of  the  surface  to  the  geometric  properties 
and  shape  of  the  surface,  and  this  relationship  makes  parametric  surfaces  well- 
suited  to  interactive  design.  Also,  many  useful  techniques,  including  subdivision 
and  local  shape  control,  have  been  developed  and  extensively  applied  in  surface 
interrogations;  see,  e.g.,  [26].  The  parametric  representation  can  be  generalized  to 
define  a  2-surface  in  Rn,  for  n  >  1. 

1.1.2  Implicit  Surfaces 

An  algebraic  surface  in  R3  is  defined  as  the  zero  set  of  an  implicit  equation 

h(x.ij.z)=  0  (1.2) 

whe.  3  h  is  a  polynomial  in  x.ij.z.  As  mentioned  before,  efficient  conversion  be¬ 
tween  parametric  and  implicit  representations  would  be  of  critical  importance  in 
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geometric  computations.  Sederberg  [61,  63]  demonstrated  that  it  is  always  possi¬ 
ble  to  implicitize  a  parametric  curve  or  surface  using  elimination  techniques  such 
as  resultant  methods  [46,  47,  49,  5S].  The  implicitization  can  also  be  based  on 
Grobner  Basis  techniques  [18,  31,  33].  Both  methods  are  fairly  expensive  and 
hence  their  practical  application  is  currently  quite  limited. 

Parameterization  or  converting  an  implicit  representation  to  an  equivalent 
parametric  representation  is  not  always  possible  since  not  all  implicit  surfaces  can 
be  expressed  as  rational  parametric  surfaces.  In  fact,  only  implicit  curves  of  genus 
zero  possess  a  rational  parameterization;  see  [1,  2,  3,  4,  33,  64]. 

i.1.3  Implicit  Surfaces  in  High-Dimensional  Space 

While  many  surfaces  can  be  formulated  quite  easily  in  three  dimensional  space, 
as  parametric  or  implicit  surfaces,  certain  surfaces  including  offsets,  Voronoi  sur- 
faces( equidistance  surfaces),  and  spherical  bk  >  can  not.  Due  to  the  possibility 
of  high  algebraic  degrees,  many  geometric  operations  on  such  surfaces  can  lead 
to  high  computational  complexity  and  numerical  instability,  and  may  require  ex¬ 
pensive  symbolic  manipulations.  As  an  alternative,  Hoffmann  proposed  in  [30.  31] 
a  surface  representation  that  is  defined  in  a  higher  dimensional  space,  with  more 
variables  but  simpler  equations.  With  this  representation,  complex  symbolic  com¬ 
putations  and  numerically  delicate  operations  can  often  be  avoided,  and  hence 
practical  implementations  can  be  realized. 

In  the  following,  we  give  the  formulation  of  offset  and  Voronoi  surfaces  in  high 
dimensional  space. 

Offset  Surfaces  As  described  in  [30],  the  r-ojfset  of  the  surface  g(x,y,z)  =  0  is 
the  set  of  points 

Offset(g, r)  =  {  p  |  d (g,p)  =  r} 

where  dig,  p)  is  the  Euclidean  distance  of  the  point  p  from  the  surface  g  =  0.  The 
offset  surface  in  general  has  two  sheets  in  real  affine  space  and  can  be  defined 
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mathematically  applying  the  envelope  theorem  as  follows: 

(x  —  u)2  +  (y  —  v)2  +  (z  —  w)2  —  r2  =  0 

g(u,v,w)  =  0 
(x  —  u,  y  —  v,  z  —  w)  •  tj  =  0 
(x  —  u,  y  —  v,  z  —  w)  •  t2  =  0 


(1.3) 


where  (u,  v,  w )  is  the  footpoint,  and  ti  and  t2  are  two  linearly  independent  tangent 
directions  of  g  —  0  to  which  the  direction  vector  (x  —  u,y  —  v,  z  —  w)  must  be 
perpendicular.  For  example,  with  =  (gw,0,  —gu)  and  t2  =  (O,#*,,—  gv),  we 


obtain 

gw(x  -  u)  -  gu(z  -  w)  =  0 
9w(V  ~  v)  ~  9v{z  -  w)  =  0 

Its  closed-form  can  be  obtained  in  principle  by  eliminating  u,u,  and  w  using  re¬ 
sultant  techniques  [13,  46,  61]  or  Grobner  basis  methods  [18,  19].  When  gw  =  0. 
(<7«m0,—  gu)  and  (0,5^,  —  gv)  become  linearly  dependent  and  hence  there  will  be 
an  extraneous  factor  representing  spheres  on  g  —  0  H  gw  =  0  in  its  closed-form 
representation.  To  eliminate  this  extraneous  factor,  we  may  adjoin 


Su(y  -v)-  gv(x  -  u)  =  0 


to  system  (1.3).  Moreover,  the  presence  of  singularities  also  introduces  additional 
extraneous  factors.  Thus,  the  formulation  in  (1.3)  is  not  faithful  in  the  sense  that 
the  natural  projection  of  the  solution  set  of  the  system  contains  points  that  are 
not  on  the  offset  surface:  see  [32].  The  offset  formulation  (1.3)  is  a  system  of  four 
equations  with  degrees  as  high  as  g's  while  its  closed-form  in  (x,  y,  x)-space  has  a 
much  higher  degree.  For  example,  when  g  is  a  quadric,  the  closed-form  of  its  offset 
may  have  degree  8.  The  computation  of  the  intersection  of  the  offset  with  other 
surfaces  is  reported  in  [31,  30].  Below,  we  define  the  offset  of  an  ellipsoid  as  an 
example. 


Example  1.1  Consider  an  ellipsoid  g{x.y,z)  =  4x2  +  y2  +  z2  —  4  =  0  which  is 
centered  at  the  origin  and  has  semiaxes  1,2.  and  2.  respectively.  The  offset  of 
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g  =  0  by  the  distance  2  is  formulated  as 

(x  —  u)2  +  {y  —  t>)2  +  (z  —  w)2  -  4  =  0 
4u2  +  v2  4-  iv2  —  4  =  0 

(1.4) 

w(x  —  u)  —  4u(z  —  w)  =  0 
w(y  —  v)  —  v(z  —  w)  =  0 

When  w  =  0,  (u;,0,  — 4u)  and  (0,u>,—  v)  become  linearly  dependent.  Hence  the 
closed-form  of  the  offset  of  g  in  (x,y,  z)-space  reveals  an  extraneous  factor  z  and 
has  degree  9.  g 


Voronoi  Surfaces  The  Voronoi  surface  of  two  given  surfaces  g(x,y,z)  =  0  and 
h(x,y,z)  =  0,  denoted  as  Vor(g,  h),  is  the  locus  of  points  equidistant  from  g  and 
h.  It  is  formally  defined  by 

Vor(<7,  M  =  {  P  €  R3  |  d(g,p)  =  d(/i,p)} 


The  Voronoi  surface  proves  useful  in  defining  constant-radius  blends  [32],  variable- 
radius  blends  [21,  30],  and  skeletons  (medial-axis  surfaces)  of  an  object  [25].  As 
shown  in  [30],  we  define  the  Vor (g,h)  as  the  common  points  of  the  offsets  from 
both  g  and  h  by  an  identical  but  unspecified  distance  r.  Thus.  Vor (g,h)  can  be 
formulated  as  eight  equations  in  ten  variables 


(x  -  u)2  4-  (y  -  u)2  +  (z  -  w)2  -  r2  =  0 

g(u,v,w)  =  0 
(x  —  u,  y  —  u,  z  —  w)  •  tt  =  0 

(x  —  u,y  —  v,z  —  w)  ■  t?  =  0 

(x  -  u)2  +  (y  -  u)2  +  (c  -  tx>)2  -  r2  =  0 

g\u,v,  w)  =  0 

(x  —  u,y  —  D.  z  —  w)  •  tj  =  0 

(x  —  u.  y  —  v.  z  —  w)  ■  t2  =  0 


(1.5) 


where  (u.v.w)  and  (u.h.  w)  are  footpoints  on  g  and  h  respectively,  and  (ti,tj) 
and  ( tx .  1 1 )  are  two  linearly  independent  tangent  directions  to  g  and  h  respectively. 


I 


By  eliminating  the  variables  u ,  v,  w,  u .  v ,  u),  and  r,  the  closed-form  of  the  Voronoi 
surface  Vor(y,  h)  can  be  obtained  in  principle. 

The  closed-form  of  Vor(<7,  h)  consists  of  two  components,  one  is  called  the  even 
Voronoi  surface  of  g  and  h  and  the  other  is  the  odd  Voronoi  surface  of  g  and  h.  The 
even  Voronoi  surface  is  the  points  that  are  either  inside  both  g  and  h  or  outside 
both  of  them.  The  odd  Voronoi  surface  consists  of  points  that  are  inside  one  and 
outside  the  other  basis  surface.  Both  components  arise  when  we  cannot  distinguish 
positive  and  negative  off ests.  Notice  that  the  elimination  methods  used  to  compute 
the  closed-form  of  (1.5)  cannot  distinguish  between  even  and  odd  Voronoi  surfaces 
and  hence  the  closed-form  must  be  the  union  of  two  surfaces.  Since  Vor(g,  h)  is 
defined  based  on  offset  formulations,  its  closed-form  may  reveal  extraneous  factors 
which  may  be  eliminated  by  adding  additional  equations:  see  [32],  To  illustrate 
the  above  formulation,  we  give  the  following  example. 

Example  1.2  Let  g  =  x2  +  (z  —  2)2  —  1  and  h  =  y2  +  (2  +  2)2  —  1  be  two  cylinders  of 
unit  radius.  Let  r  be  the  common  radius  of  the  spheres  centering  on  g  and  h,  and 
(j,y,z)  be  a  point  equidistant  from  g  and  h  at  (u,v,u;)  and  (u,u,  u>)  respectively. 
At  ( u,v,w )  on  g  =  0,  t[  =  (0,1,0)  and  t2  =  (-(u;  -  2),0,u)  are  two  linearly 
independent  tangent  vectors.  At  (u.v.w)  on  h  =  0,  two  linearly  independent 
tangents  are  t!  =  (1,0,0)  and  t2  =  (0.  — ( tZ>  4-  2),  v).  We  then  have  the  following 
system  of  eight  equations  representing  Vor (g,h) 

(x  —  u)2  +  (ij  —  v)2  +  (z  -  w)2  —  r2  =  0 
u2  4-  (w  —  2)2  —  1  =  0 

u{z  —  w)  +  (2  —  tr)(x  -  u)  =  0 

y  -  v  =  0 

(1.6) 

(x  —  u)2  +  (f/  —  'M2  *  —  w)2  —  r2  =  0 

■2  ~  i  ff  -  2  j 2  -  1  =  0 
f’(  -  —  'Cl  -  ■  v  —  2  —  v)  =  o 


.f 


0 


s 


After  eliminating  u,  v ,  w ,  u,  u,  u>,  r  from  (1.6),  we  obtain 

(x2  —  y2  —  8z)(48z2  4-  16y2z  —  16x2r  +  y4  —  2i2y2  —  8y2  +  x4  —  8x2  —  48) 

where  the  first  factor  is  the  even  Voronoi  surface  of  g  and  h  and  the  second  one 
is  the  odd  Voronoi  surface.  Note  that  the  extraneous  factors  do  not  appear  in  the 
closed-form  of  Vor(g,h)  due  to  the  existence  of  two  linearly  independent  tangents 
at  every  point  of  the  cylinders.  ■ 

Blending  Surfaces  Given  two  surfaces  g(x,  y,  z)  =  0  and  h{x,  y,  z)  =  0.  a  blend¬ 
ing  surface  is  a  surface  that  intersects  both  surfaces  tangentially  along  two  curves. 
A  constant-radius  blend  is  a  blending  surface  that  has  circular  cross-sections  of 
fixed  radius.  A  variable-radius  blend  is  a  blend  whose  circular  cross-sections  are 
of  variable  radius.  A  constant-radius  blend  of  g  and  h  can  be  formulated  as  the 
envelope  of  the  family  of  spheres  of  constant  radius  r  whose  centers  are  constrained 
to  lie  on  Offsetfy,  r)  fl  Offset(A,  r).  Formalized  in  algebraic  terms,  it  is  the  zero  set 
of  ten  equations  in  twelve  variables;  see  [30].  A  variable-radius  blend  of  g  and  h 
is  the  envelope  of  the  family  of  spheres  that  have  centers  lying  on  Vor (g,h)  fl  p, 

where  p  is  a  reference  surface,  and  have  radii  such  that  each  sphere  touches  both 

g  and  h:  see  also  [21,  30]. 

When  the  basis  surface  of  the  offset  is  in  parametric  form 
x  =  x(u,v),  y  =  y(u,v\,  z  =  z(u,v) 

the  algebraic  formulation  results  in  3  equations  in  5  variables  as  follows 

(x  -  x(u,u)j2  +  (y  -  y(u.t-))2 (z  -  ;(u.u))2  -  r2  =  0 

(x  x( u,  v),  y  — *  y(u,  v),  z  —  z(u,  v) )  ■  (xu(u,v),yu(u,v),zu(u,u))  =  0  G-") 

(x  -  x(u.v),y  -  y(u,v),z  -  z(u.y))  •  fx,.(  u.  v).  yv(  u.  v ).c„(u.  v))  =  0 

where  the  subscripts  denote  partial  differentiation.  For  Voronoi  surfaces  and  blends 
involving  parametric  basis  surfaces,  the  formulation  is  closely  analogous  to  offsets. 
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This  algebraic  formulation  of  offsets,  Voronoi  surfaces  and  blends  results  in  2- 
surfaces  in  Rn.  n  >  3,  which  are  generally  defined  as  the  zero  set  of  the  following 
system  of  polynomial  equations: 


/i(xi,x2,...,xn)  =  0 

f2(Xl,X2,....Xn)  =  0 


(1.8) 


f m  ( ^1  *  %2  i  •  •  •  *  Xji  )  —  0 

where  /,  €  A'fij,  x2, . . .  ,x„]  and  m  (>  n  -  2)  is  normally  n  -  2.  We  denote 
system  (1.8)  in  matrix  form  as  F(x)  =  0  and  the  zero  set  of  F(x)  =  0  as  Sp. 
Recall  that  extra  equations  are  sometimes  required  in  order  to  eliminate  extraneous 
factors  introduced  by  linearly  dependent  tangent  vectors  in  the  offset  formulation 
as  shown,  e.g.,  in  (1.3).  In  the  rest  of  the  thesis,  we  assume  m  =  n  —  2.  For  the 
cases  of  m  >  n  —  2  the  modifications  needed  in  the  proposed  computations  are 
routine.  It  is  worth  remarking  that  the  algebraic  formulation  of  offsets,  blends, 
and  Voronoi  surfaces  given  here  all  have  the  property  that  a  2-surface  is  defined 
in  Rn,  where  n  >  3,  but  its  projection  into  a  certain  subspace  is  wanted.  In  this 
thesis,  we  assume  that  the  (ii,  12,  i3)-space  is  this  subspace. 

In  principle,  from  system  (l.S),  the  last  n  —  3  variables  of  x  can  be  eliminated. 
This  computation  reduces  F(x)  =  0  to  a  single  equation 


/(zj,z2,z3)  =  0  (1.9) 

in  (xi,  Xj,  z3) -space  with  its  zero  set  denoted  as  S/.  However,  the  elimination 
process  is  not  practical.  Hence  the  closed-form  representation  of  F(x)  =  0  in 
(zt,  Z2,  z3)-space,  f(x i,Z2,x3)  =  0,  is  often  unobtainable  in  practice.  In  general 
S/  is  the  natural  projection  of  Sp.  Note  that  Sf  might  contain  more  points  than 
the  natural  projection  of  Sp  since  the  projection  of  Sp  need  not  be  an  algebraic 
variety:  see  [69]. 
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1.2  Implicit  Approximations  of  Parametric  Curves  and  Surfaces 

A  recurring  operation  in  solid  modeling  is  the  evaluation  of  surface  intersections 
[55].  If  both  surfaces  are  given  parametrically,  the  two  major  approaches  given  the 
greatest  prominence  in  the  literature  are  subdivision  and  substitution  methods. 

In  the  subdivision  method,  e.g.,  [3S.  40,  41.  42,  53],  both  surfaces  are  recursively 
subdivided  in  the  vicinity  of  their  intersection.  The  subdivision  results  in  an 
adaptive  piecewise  linear  approximation  of  both  surfaces  and  their  intersection. 
Among  the  advantages  of  the  method  we  mention  its  robustness  and  its  potential 
for  locating  all  intersection  branches.  A  major  drawback  of  the  subdivision  method 
is  the  large  volume  of  data  it  creates,  which  slows  it  down  in  areas  of  high  surface 
curvature. 

In  the  substitution  method,  e.g,  [27,  44,  59.  67,  68],  one  of  the  surfaces,  Si,  is 
converted  to  implicit  form  F,  and  the  parametric  form  of  S2  is  substituted  into  F 
resulting  in  an  implicit  algebraic  curve  /  in  the  parameter  space  of  S2.  This  curve 
/  is  in  birational  correspondence  with  the  intersection  of  S\  and  S2  in  xyz-space, 
and  thus  serves  as  an  accurate  representation  of  the  intersection.  Major  difficul¬ 
ties  of  the  substitution  method  limit  its  utility  in  practice.  There  are  two  general 
methods  for  implicitizing  a  parametric  surface.  The  first  method  is  based  on  elim¬ 
ination  theory  [61]  and  does  resultant  computations.  It  is  expensive  and  generates 
extraneous  factors  whose  detection  is  a  delicate  problem,  see  also  Section  2.3.  The 
second  method  for  implicitization  is  based  on  Grobner  basis  techniques  [IS].  It  is 
’.so  fairly  expensive  and  requires,  moreover,  rational  coefficients  in  the  description 
of  S\.  Another  difficulty  with  the  substitution  method,  less  prominently  pointed 
out  but  well-known  [51],  is  that  the  substitution  itself  can  be  numerically  unstable, 
and  is  a  nontriviai  algorithmic  task  when  desiring  efficiency  and  accuracy.  Some 
authors  have  suggested  the  use  of  rational  arithmetic  for  this  reason  [2S],  thus 
further  adding  to  the  computational  load  of  the  approach. 
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In  Chapter  2,  we  provide  a  middle  ground  by  deriving  a  local  implicit  approx¬ 
imation  of  rational  or  polynomial  parametric  curves  and  surfaces  with  low-degree 
implicit  forms.  In  the  context  of  subdivision  techniques,  such  approximations  have 
the  potential  of  reducing  the  number  of  generated  surface  approximants,  because 
we  are  not  restricted  to  linear  approximants  only.  In  the  context  of  substitution 
methods,  the  approximations  avoid  the  high  cost  of  implicitizing  a  parametric 
curve  or  surface,  and  provide,  moreover,  irreducible  approximants. 

Since  the  distribution  of  a  preliminary  version  of  [23],  a  number  of  related 
investigations  have  been  developing  and  applying  similar  ideas.  Bajaj  and  Ihm  [15] 
apply  a  technique  that  is  analogous  to  ours  to  the  problem  of  designing  blending 
surfaces  and  prove  results  on  minimum  degree  blends  satisfying  certain  constraints. 

Previously,  local  explicit  approximations  to  integral  parametric  curves  and  sur¬ 
faces  have  been  proposed  in  [4S].  An  approximant  of  the  form 

2  =  f{x,y)  =  or  V  =  f[x)  =  5Za«x' 

is  constructed,  for  surfaces  and  curves.  Recurrence  formulas  were  also  derived  for 
the  coefficients  of  /.  Bajaj  [11]  extends  this  method  using  power  series  compo¬ 
sition  and  inversion  techniques  together  with  rational  Pade  approximations.  In 
our  experience,  a  local  explicit  approximation  is  less  favorable  than  a  local  im¬ 
plicit  approximation.  In  fact,  while  a  quadratic  explicit  approximation  to  a  curve 
achieves  second  order  contact  at  the  point  at  which  it  is  constructed,  a  quadratic 
implicit  approximation  achieves  fourth  order  contact.  For  curves,  the  order  of 
contact  grows  linearly  with  the  degree  of  the  explicit  approximation,  whereas  the 
order  of  contact  of  the  implicit  approximation  has  a  quadratic  growth  in  the  de¬ 
gree.  Thus,  much  lower  degree  approximations  suffice.  Note,  however,  that  for  an 
implicit  approximation  of  degree  n.  O(n)  coefficients  can  be  chosen,  whereas  for  a 
degree  n  explicit  approximation  only  O(n)  coefficient  are  available. 

In  general,  local  explicit  approximation  can  only  approximate  curves  or  sur¬ 
faces  locally  no  matter  how  high  a  degree  of  approximant  is  used.  This  is  due 
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to  the  asymmetry  introduced  by  making  one  variable  an  explicit  function  of  the 
other(s).  For  instance,  a  circle  cannot  be  completely  approximated  with  a  single 
explicit  approximant.  In  contrast,  our  approximants  are  capable  of  approximating 
curves  or  surfaces  not  only  locally  but  also  globally  in  the  sense  that  the  radius  of 
convergence  increases  when  the  degree  of  approximation  increases,  and  the  exact 
implicitization  can  be  finally  derived  when  the  degree  of  approximation  is  equal  to 
the  degree  of  the  given  parametric  curve  or  surface. 

For  a  properly  parameterized,  rational  curve  r(t),  e.g.,  [60,  62],  of  degree  m 
containing  the  origin,  we  seek  an  implicit  curve  g{x,y )  =  0  of  degree  n  <  m  that 
approximates  r(£)  at  the  origin.  The  idea  is  to  set  up  the  polynomial  g{x,y)  with 
symbolic  coefficients  e,j  and  substitute  r (t)  into  g(x,y)  with  result 

nm 

5(r(0)  = 

t=i 

where  the  a,  are  linear  combinations  of  the  e,}.  We  require  that  a  certain  number 
of  the  q,  vanish.  With  one  of  the  coefficients  being  set  to  one  and 

e*i  =  0,  a2  =  0,  ...,  a,  =  0 

for  some  s  the  e,j  are  determined  and  an  implicit  approximation  is  obtained  that 
has  contact  of  order  s  with  r(f)  at  the  origin. 

We  have  derived  a  recurrence  for  computing  the  a,  directly  from  r (t)  without 
explicit  substitutions  and  shown  that,  when  n  <  m, 

•  the  coefficient  matrix  of  an,  Q2, . . . ,  anm  has  rank  that  is  equal  to  the  number 
of  unknown  coefficients  in  g(x,y). 

•  the  degree  n  local  implicit  approximation  g(x,y)  of  r(t)  at  the  origin  is  irre¬ 
ducible  when  the  origin  is  a  regular  curve  point. 

The  method  derived  here  has  the  following  merits  compared  to  local  explicit 
approximations  proposed  in  (48.  11]: 
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•  It  works  for  polynomially  and  rationally  parameterized  curves  and  surfaces. 

•  It  yields  meaningful  results  for  many  types  of  singularity. 

•  The  order  of  contact  grows  quadratically  for  curves  and  faster  than  linear  for 
surfaces. 

•  The  implicit  approximants  approximate  curves  and  surfaces  not  only  locally 
but  also  globally. 

The  algorithm  for  computing  the  implicit  approximation  of  a  properly  parameter¬ 
ized  rational  surface  is  basically  similar  to  the  one  for  the  curve  case  except  that 
in  the  surface  case  additional  safeguards  are  incorporated  to  ensure  irreducibility. 

1.3  Local  Approximations  of  2-D  Surfaces 

As  described  in  Section  1.1.3,  certain  surfaces,  including  offset,  blending  and 
Voronoi  surfaces,  cannot  be  easily  defined  in  the  conventional  3-D  space.  In  con¬ 
trast,  these  surfaces  can  be  formulated  mathematically,  with  the  theory  of  en¬ 
velopes,  in  higher  dimensional  space  in  a  straightforward  manner.  As  the  result, 
such  surfaces  are  generally  2-dimensional  surfaces  in  Rn,  n  >  3,  and  are  defined 
by  (l.S)  or  in  a  matrix  form  F(x)  =  0.  Although  the  exact  closed-form  representa¬ 
tion  of  such  a  surface  /(xt,  x2l  x3)  =  0  could  be  derived  in  principle  by  elimination 
methods  such  as  Gobner  basis  [18,  19]  or  resultant  techniques  [46,  61,  13],  it  is 
often  not  feasible  to  do  so  due  to  the  high  complexity  of  these  methods.  Thus, 
as  mentioned  in  [32],  the  viable  approaches  to  interrogating  such  surfaces  seem  to 
be  (a)  approximate  the  surface  locally  and  interrogate  the  approximant.  or  (b)  in¬ 
terrogate  the  higher-dimensional  representation.  In  Chapter  3,  we  have  developed 
computational  schemes  for  the  local  geometry  of  f(xi ,  x2,  x3)  =  0  and  have  investi¬ 
gated  techniques  that  construct  implicit  approximants.  explicit  approximants.  and 
parametric  approximants  to  the  surface  in  (xi ,  x2,  X3)-space  described  by  a  system 
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of  algebraic  equations  (1.8).  These  approximants  have  the  forms  x3  =  y>(x!.x2), 
and  (xx  =  0i(s,f),x2  =  <p2(s,t),x3  =  <z>3(s,f)),  respectively. 

Since  /  is  unknown  and  its  derivation  is  often  impractical,  methods  of  comput¬ 
ing  the  normal,  tangent,  and  normal  curvatures  of  /  =  0  at  a  surface  point  in  R3 
from  the  information  depending  only  on  system  (1.8)  have  been  developed.  More¬ 
over,  using  a  variant  of  the  Three  Tangent  Theorem  given  by  Pegna  and  Wolter 
[52],  a  method  has  been  developed  for  computing  a  degree  two  local  implicit  ap- 
proximant  of  /  =  0. 

The  implicit  function  theorem  ensures  that  system  ( l.S)  determines  m  =  n  —  2 
components  of  x  =  (xt,...,xn)  as  functions  of  the  remaining  2  components  in  a 
neighborhood  of  any  surface  point  x°  on  which  the  differential  D F(x°)  has  rank 
n  —  2.  The  unknown  coefficients  of  explicit  functions  can  be  calculated  by  means 
of  the  chain  rule  and  a  linear  system  solver;  however,  it  is  algebraically  tedious  to 
do  so.  When  we  assume  that  x°  is  the  origin,  a  recursive  formula  has  been  derived 
which  presents  the  computation  in  a  more  convenient  manner.  Thus,  the  recursive 

formula  computes,  without  loss  of  generality,  x<  =  ^i(xi,x2),  i  =  3, _ n.  It  is 

clear  that  x3  =  t/;3(xi,x2)  is  a  local  explicit  approximation  of  /  =  0. 

For  a  given  system  (l.S)  and  a  regular  point  x°  on  it.  there  exists  a  neighbor¬ 
hood  of  x°  in  which  the  parametrically  defined  solution  x,  =  o,(u,v),i  =  1,2 . n. 

can  be  found  such  that  x°  =  (<z>i(0,0) . o„(0, 0)).  The  first  three  coordi¬ 

nate  functions  (<t>i(u,v)td>2(u,v),<p3(u,v))  so  computed  constitute  therefore  a  local 
parametrization  of  /  =  0.  To  compute  the  parametric  solution,  we  substitute 

symbolically  x,  with  <p,(it,u),  i  =  l,...,n,  in  ail  polynomials  /,, i  =  1 . n  —  2. 

Then  we  compute  the  Taylor  expansion  of  the  resulting  polynomials  in  u  and  v. 
By  requiring  that  the  coefficient  of  uJvk  is  identically  zero  for  1  <  j  +  k  <  /,  a 
series  of  linear  systems  are  obtained.  The  solutions  of  the  linear  systems  define 
a  degree  /  approximation  of  the  parametrically  defined  solution.  The  first  three 
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component  functions  xx  =  o1(u,u).X2  =  02(11.  u),x3  =  <p3(u,u)  so  obtained  com¬ 
prise  the  degree  l  parametric  approximant  of  /  =  0.  This  derivation  is  analogous 
to  the  derivation  of  approximants  of  space  curve  proposed  in  [14]. 

1.4  Piecewise  Linear  Approximations  of  2-D  Surfaces 

With  the  availability  of  piecewise  linear  approximations  of  surfaces,  certain 
computations  become  practical;  in  particular,  surface  rendering  can  take  advantage 
of  hardware  capabilities  and  we  do  not  have  to  resort  to  expensive  ray  casting 
for  visualization  purposes.  Based  on  subdivision,  the  PLA  of  major  parametric 
surfaces  in  CAGD,  including  Bernstein-Bezier  surface  and  B-spline  surfaces,  have 
been  extensively  studied  and  used  with  great  success  in  surface  interrogations 
[10,  16,  38,  41,  42].  However,  it  seems  that  much  less  research  has  been  done  in 
the  literature  to  the  PLA  of  implicitly  defined  surfaces  [17,  9,  7,  S,  56,  57,  12]. 

Bloomenthal  [17]  has  proposed  an  algorithm  for  computing  PLA  of  an  implicit 
surface  g(x,y ,  z)  =  0  based  on  space  subdivision  using  octrees.  The  algorithm 
starts  with  an  octree  that  bounds  the  surface  portion  of  interest,  and  then  decom¬ 
poses  recursively  those  octrees  that  intersects  the  surface  until  the  surface  portion 
inside  the  octree  is  sufficiently  close  to  a  plane.  The  surface  function  is  evaluated 
at  the  corners  of  an  octree  cube,  and  from  these  values  is  determined  the  point 
at  which  the  surface  intersects  an  edge  of  the  cube,  by  linear  interpolation.  The 
surface  points  on  the  edges  of  each  cube  are  ordered  to  form  a  convex  polygon 
each  of  whose  edges  lies  in  a  cube  face.  One  disadvantage  of  the  vertex  evaluation 
strategy  is  that  the  negatively  signed  corners  of  a  cube  cannot  be  separated  from 
the  positively  signed  corners  by  a  single  plane  in  all  cases.  Thus,  vertex  evaluation 
on  a  cube  may  result  in  ambiguities  that  could  produce  more  than  one  polygo¬ 
nal  approximation.  In  this  case,  further  subdivision  is  usually  required  but  does 
not  necessarily  resolve  the  ambiguity.  Due  to  the  nature  of  space  decomposition. 
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this  approach  is  capable  of  approximating  all  the  surface  components  provided  the 
initial  cube  has  been  properly  situated. 

In  [9,  7,  8.  6],  a  simplicial  continuation  algorithm  is  presented  for  obtaining  a 
PLA  to  a  component  of  an  implicitly  defined  2-surface  in  Ft".  The  algorithm  starts 
at  a  regular  surface  point  x°  and  generates  n-simplices  in  a  sequential  ordering 
which  spiral  outward  from  x°.  The  approximate  zero  points  on  the  (n  —  2)-faces 
of  each  n-simplex  are  computed  by  an  affine  map  and  can  be  further  refined  using 

Newton  iteration.  Note  that  a  n-simplex  a  =  [v0,  Vt . vn]  is  the  set  of  points 

in  Rn  that  is  the  convex  combination  of  v0,  . vn.  A  n-simplex  a  is  said  to 

be  transversal  if.  roughly  speaking,  there  exists  an  (n  —  2)-face  with  which  the 
2-surface  intersects  transversely.  The  standard  vector  labeling  of  v  6  Rn  induced 
by  F 


lF(v) 


1 

F(v) 


is  used  to  formally  define  transversality  as  f  lu/.vs.  An  (n— 2)-face  [v0,  V! . vn_2j 

is  said  to  be  completely  labeled  with  respect  to  the  vector  labeling  lp  if  the  labeling 
matrix 


A  = 


1 


[  F(v0)  F( vn_2)  J 

has  a  lexicographically  positive  inverse,  i.e..  the  first  nonzero  element  in  each  row 
of  A"1  is  positive.  An  n-simplex  a  is  transversal  if  there  exists  an  (n  -  2)-face 
r  C  <7  which  is  completely  labeled  with  respect  to  the  vector  labeling  1F. 

For  a  n-simplex,  the  affine  map  :  R"  — *  R"-2  is  a  linear  map  which 

interpolates  Fix)  on  the  vertices  of  a.  i.e..  HJv,)  =  F(v,)  for  i  =  0.1 . n. 

When  cr  is  transversal  cr  n  H7 1  ( 0 )  serves  as  an  approximation  to  the  2-surface 
inside  cr. 

Let  T  represent  the  collection  of  oil  transversal  simplices  in  a  compact  domain 
and  the  piecewise  linear  approximation  Ht  of  F(x)  relative  to  T  be  defined  as  H„ 
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when  restricted  to  <7  for  every  a  €  T.  Then, 

Ht  (0)  =  {  HJ^O)  |  <7  6  T} 

is  the  PLA  of  the  surface  component  that  contains  x°.  Newton  iteration  can  be 
used  to  refine  points  on  H^l(0)  to  surface  points.  An  error  bound  for  ||F(x)||, 
where  x  6  Hi^O),  in  terms  of  the  mesh  size  is  provided  and  it  is  shown  that 
H^x(0)  approximates  F-1(0)  quadratically  in  the  mesh  size  [5]. 

While  Bloomenthal’s  method  is  more  suited  to  implicit  surfaces  in  R3,  the 
simplicial  continuation  method  is  designed  for  general  2-surfaces  in  Rn.  However, 
the  direct  application  of  these  methods  to  computing  PLA  of  offsets,  Voronoi 
surfaces,  and  blends  will  have  a  space  and  time  complexity  that  is  exponential  in 
the  dimensionality  of  the  space. 

Rheinboldt  [56,  57]  has  presented  a  continuation  method  that  maps  a  reference 
triangulation  £  on  Rp  to  a  p- manifold  M  in  Rn,  n  >  p  >  1,  and  hence  produces  a 
triangulation  on  M.  The  method  consists  of  two  major  computation  schemes.  A 
moving  frame  algorithm  is  given  to  derive  orthonormal  bases,  i.e.  local  coordinate 
systems  of  the  tangent  space  TZ{M ),  that  vary  continuously  with  their  point  of 
contact  x  on  A/.  A  triangulation  algorithm  uses  the  orthonormal  bases  produced 
bv  the  moving  frame  algorithm  to  map  the  vertices  of  the  reference  triangulation 
onto  the  tangent  space  x  +  TZ(M)  corresponding  to  appropriate  points  x  on  M. 
Thereafter,  Gauss-Newton  iteration  is  applied  to  project  these  triangulations  from 
x  +  Tr(  A/)  onto  M.  At  each  step  of  the  triangulation  algorithm,  a  reference  vertex 
c  of  the  reference  triangulation  £  and  the  corresponding  point  x  6  A/  are  selected. 
The  center  £  is  associated  with  a  set  T(£)  of  vertices  of  £  that  can  be  mapped 
onto  A/.  Then  the  reference  vertices  in  )  that  have  not  been  processed  are 
mapped  onto  x +  TX(A/)  and  projected  onto  M  by  the  Gauss-Newton  iteration.  To 
proceed,  a  reference  vertex  in  T(^)  that  has  been  processed  is  selected  as  the  next 
c  and  the  same  computation  is  applied.  The  points  computed  on  M  inherit  the 
connectivity  structure  of  £  which  m  turn  ;;:<iuces  a  triangulation  on  A/.  Note  that 
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the  triangulations  induced  by  two  overlapped  T(£i)  and  T(£2)  are  made  compatible 
by  demanding  that  the  reference  vertices  in  T(fi)  fl  T(£2)  are  processed  only  once. 

In  Chapter  4,  we  propose  an  algorithm  that  is  capable  of  computing  PL  A  for 
2-surfaces  projected  into  R3  but  defined  algebraically  in  Rn,  n  >  3. 

1.5  Some  Preliminaries 

A  rational  plane  curve  r(t)  can  be  given  as  the  pair  (x(t),  y(f)),  where  x(t)  and 
y{t)  are  rational  functions  of  t.  The  curve  points  are  all  points  (x(f).y(<))  on  the 
plane.  The  curve  is  properly  parameterized  if  for  all  but  finitely  many  curve  points 
p  we  have  p  =  (x(f),y(t))  with  a  unique  value  of  t.  When  a  parametric  curve  is 
not  properly  parameterized,  there  exists  a  rational  nonlinear  function  s(t)  such 
that  x(t)  =  x’(s(t))  and  y(t)  =  ym(s(t))  for  some  rational  functions  x*  and  y*.  We 
assume  in  this  thesis  that  all  parametric  curves  are  properly  parameterized  and 
note  that  a  parametric  curve  is  always  irreducible.  For  methods  to  detect  improper 
parameterization,  see  Sederberg  [62]. 

The  degree  of  a  rational  parametric  curve  is  the  highest  degree  of  the  numerator 
or  the  denominator  polynomial,  assuming  both  x(t)  and  y(t)  have  been  written 
with  a  common  denominator.  The  implicit  equation  /(x,y)  of  the  rational  curve 
r(t)  is  a  lowest  degree  polynomial  in  x  and  y  satisfying  /(x(t),yU))  =  0.  It  is 
unique  ud  to  a  multiplicative  constant.  If  r(t)  has  degree  m.  then  so  does  f(x.y); 
see,  e.g.,  [46,  49]. 

As  with  parametric  curves,  a  parametric  surface 

P (s.t)  =  (x(s,  t),  y(«,  t)yz(s.  t)) 

can  be  improperly  parameterized  if  there  are  nonlinear  rational  functions  u(s.t) 
and  t)  such  that 


P  (s.t)  =  (x‘(u(s.t),v($,t)),  ym{  ul. s.t ).  v(s.t)),  :m{u(s.t)*v(S't))) 


for  some  x",y *  and  z“.  In  that  case,  there  is  a  many-to-one  correspondence  be¬ 
tween  the  parameter  values  and  the  surface  points.  P(s,f)  is  properly  parame¬ 
terized  if  this  correspondence  is  one-to-one  except,  possibly,  on  a  one  dimensional 
set  of  points.  We  also  assume  that  all  parametric  surfaces  are  properly  parame¬ 
terized.  For  a  parametric  surface  described  by  rational  functions  of  total  degrees 
m,  there  always  exists  an  irreducible  implicit  equation  f(x,y,z)  =  0  satisfying 
/(x(s,  t),  y(s,  t),  z(s,  t))  =  0,  and  /  is  unique  within  a  constant  factor,  see,  e.g., 
[22].  Moreover,  /  has  degree  at  most  m2. 

A  polynomial  of  degree  k  in  the  variables  xj,x2,...,xn  is  denoted  fk( xlf ...,  xn) 
whenever  we  wish  to  stress  the  degree.  The  gradient  of  /  =  0  at  the  point 
x  =  (i1,i2,...,in)  is  the  vector  V/  =  (/x,  ,/X2, . . . , /Xn),  where  the  partials  are 
evaluated  at  x.  The  Hessian  of  /  is  the  symmetric  matrix 


where  subscripts  denote  partial  differential. 

A  point  x°  =  (x°,  •••,x°)  is  said  to  be  non-singular  or  regular  on  /  if  the 
gradient  of  /  at  x°  is  not  null:  otherwise  the  point  is  singular. 

Given  a  system  of  algebraic  equations  (1.8),  the  m  x  n  matrix 


DF(x°) 


is  called  the  differential  of  F  at  x°  =  (x°,  •  •  •  ,  x°).  A  point  x°  is  said  to  be  non¬ 
singular  or  regular  on  the  2-surface  F  if  D F(x°)  has  rank  n  —  2.  That  is.  the 
gradients  V/i, . . . ,  V/m  form  a  normal  space  of  dimension  n  —  2.  When  m  =  n  —  2, 

this  is  equivalent  to  saying  that  x°  is  non-singular  on  all  /,,  i  =  1,2 . n  —  2. 

and  the  gradients  V/x, ....  V/m  are  linearly  independent.  x°  is  singular  if  it  is  not 
a  non-singular  point. 
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1.6  Thesis  Organization 

In  this  first  chapter,  various  types  of  surface  encountered  in  CAGD  and  solid 
modeling,  including  parametric  surfaces,  implicit  surfaces,  and  2-surfaces  in  high¬ 
dimensional  space,  have  been  presented,  and  the  problems  of  approximating  curves 
and  surfaces  and  their  applications  to  surface  interrogations  have  been  discussed. 
Moreover,  several  local  and  global  approximation  techniques  that  will  be  proposed 
in  this  thesis  have  been  sketched. 

In  Chapter  2,  a  method  for  computing  implicit  approximations  to  parametric 
curves  or  surfaces  is  proposed  that  might  circumvent  the  difficulty  of  globally  im- 
plicitizing  parametric  curves  or  surfaces.  Chapter  3  deals  with  the  local  geometry 
and  local  approximations  of  2-surface  in  high-dimensional  space.  In  Chapter  4,  an 
algorithm  is  developed  for  computing  the  piecewise  linear  approximation  of  offsets, 
Voronoi  surfaces,  and  blends.  Chapter  5  provides  a  summary  of  this  research  and 
comments  on  some  future  research  directions. 
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2.  IMPLICIT  APPROXIMATION  OF  CURVES  AND  SURFACES 

A.  method  is  described  for  finding  an  implicit  approximant  to  a  parametric 
curve  or  surface  in  a  neighborhood  of  a  point.  The  method  works  for  both  poly- 
nomially  and  rationally  parameterized  curves  and  surfaces,  and  achieves  an  order 
of  contact  that  can  be  prescribed.  In  the  case  of  nonsingular  curve  points,  the  ap¬ 
proximant  must  be  irreducible,  but  in  the  surface  case  additional  safeguards  have 
been  incorporated  into  the  algorithm  to  ensure  irreducibility.  The  method  also 
yields  meaningful  results  for  many  types  of  singularity.  The  chapter  is  organized 
as  follows.  In  Section  2.1,  we  describe  the  method  for  polvnomially  and  rationally 
parameterized  curves.  Section  2.2  presents  the  surface  case.  In  Section  2.3.  we 
comment  briefly  on  some  theoretical  connections  between  the  method  we  propose 
here  and  several  resultant  formulations  found  in  the  literature. 

2.1  Local  Implicit  Approximation  of  Parametric  Plane  Curves 

We  seek  an  implicit  curve  g{x,y)  =  0  that  approximates  the  parametric  curve 
r(t)  =  (x(t),  y{t))  at  the  origin  subject  to  a  prescribed  order  of  contact.  The  idea  is 
to  set  up  a  polynomial  g(x,y)  of  sufficiently  high  degree  with  symbolic  coefficients 
et;.  Then,  a  system  of  linear  equations  with  unknowns  e,;  is  formulated  and  solved. 

The  linear  system  is  obtained  by  substituting  r(£)  into  g(x,y).  The  result  is 

g{x{t),y(t))  =  atktk 

where  the  ak  are  linear  combinations  of  the  e,r  We  require  that  a  certain  number 
of  the  ak  vanish.  With 


O  1  =  0.  02=0 . Qa  =  0 


for  some  s,  an  implicit  approximation  is  obtained  that  has  contact  of  order  s  with 
r(t)  at  the  origin.  The  approach  depends  on  the  following  details: 

1.  There  is  a  recurrence  for  deriving  the  linear  system  directly  from  r(t)  without 
explicit  substitutions.  This  recurrence  is  derived  in  Section  2.1.1. 

2.  Assume  that  the  degree  n  of  the  approximation  is  smaller  than  m,  the  degree 
of  r(t).  There  is  a  function  <^(n)  that  determines  the  order  of  contact  that 
g{x,y)  can  achieve.  This  function  is  obtained  by  analyzing  the  rank  of  the 
linear  system  in  Section  2.1.2. 


In  Section  2.1.3  we  discuss  the  error  behavior  of  the  implicit  approximation,  and. 
in  Section  2.1.4,  we  present  several  experiments. 

Let 


r(t)  =  (x(t),y(t))  =  ~r) 

w(t)  w{t) 


be  a  properly  parameterized  rational  curve  of  degree  m  containing  the  origin,  where 


p(o  =  £a«^’  9(0  =  w(o  =  zicit* 

1=1  1=1  1=0 

We  assume  that  am  and  bm  are  not  both  zero,  and  that  c0  ^  0.  From  [61,  63],  we 

know  that  there  exists  an  irreducible  polynomial  fm(x.y)  =  0  of  degree  m  such 

that 


fm(x(t),y(t))  =  0 


Let  gn{x,y)  =  e,jX'yJ  =  0  be  a  degree  n  implicit  curve  containing  the 

origin.  Since  gn( x,y)  =  0  and  ~/gn{x,y)  =  0,  where  7  #  0.  are  the  same  curve. 
g'ix.  y>  =  0  has  ^(n)  =  (n2  4-  3n  —  2)/2  coefficients  on  which  the  curve  depends. 
Let  Gn(j.!/,c)  be  the  homogeneous  form  of  gn(x,y).  Substitution  yields 

n,  p(t)  <?(0  _  Gn(p{t).  q(t ),  u>(  t))  _ 

9  wit)'  w(ty  ( w(t))n  ~  ( w{t))n 


where  the  a,  are  linear  combinations  of  t he  e,;. 
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We  look  for  an  implicit  form  yn(x,y)  =  0  of  degree  n  <  m  approximating  r(t) 
at  the  origin,  and  the  method  of  deriving  should  work  whether  w{t)  =  1  or  not.  i.e.. 
irrespective  whether  the  curve  is  parameterized  polynomially  or  rationally.  The 
following  simple  example  demonstrates  the  approach. 

Example  2.1  Consider  cQ(t)  =  (x(t),y{t))  =  (p(t)/w(t),q(t)/w(t)),  where  p(t)  = 
2 t3  +  t2  —  3 t,  q{t)  =  t3  —  t2  —  2 1,  and  w(t)  =  t2  4-  4t  +  5.  c0(t)  is  a  properly 
parameterized  plane  curve  containing  the  origin.  Lety2(x,y)  =  eiox+eoiy  +  e20x2-|- 
en xy  +  e02 y2  be  a  degree  2  curve  containing  the  origin  with  symbolic  coefficients. 
Substituting  x{t)  and  y{t)  into  g2  yields  g2(x{t),y{t))  =  (IlLi  where 

=  — 15eio  —  lOeoi 

ct  2  =  —  i  Cio  —  I3e0i  +  9e20  +  6en  +  4eo2 
0:3  =  lleto  —  eoi  —  6e2o  +  en  +  4e02 
<*4  =  9eio  +  3e0i  —  lle2o  —  Sen  —  3eo2 
Q5  =  2e10  +  eoi  +  4e2o  —  eu  —  2eo2 
ate  =  4e2o  +  2eu  +  C02 

By  requiring  e10  —  1  =  O.ai  =  0,a2  =  0,a3  =  0,  and  q4  =  0,  we  can  solve  for 
the  unknown  coefficients  e,_,.  The  resulting  g2  approximates  Co  at  the  origin  with 
contact  of  order  four.  I 

2.1.1  A  Recurrence  for  ak 

Since  gn(  x,  y)  =  gn~l(x,  y)  4-  H,+;=n  eqx'y-7,  the  homogeneous  form  of  gn(x.y) 
can  be  written  as 


G"(x.y,c)  =  :Gn  *(x.y,c)+  £  e.jxV  (2-1) 

In  the  following,  let  a£-1  and  aq  denote  the  coefficient  of  tk  in  Gn_1(p(t),  q{t).  tr(  t)) 
and  Gn( p{t).q(t),  w(t)),  respectively.  It  is  clear  that  a£  can  be  derived  from  the 
a,1'1,  i  =  1.2 . k,  because  of  (2.1 ). 
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We  define  (a(i))i  and  (&(j));  as  in  [48],  setting 

m  tm 

(p(0)‘  =  (52  =  H(a(*))<*' 

i=i  i=i 


jm 


and  similarly, 

my  =  (f><')'  =  B‘u>v' 

<=1  *=j 

A  recurrence  to  compute  the  (a(i));  and  (b(j))/  is  derived  as  follows,  see  also  [48]. 
By  definition  of  (a(i  4-  l))j,  it  follows  that 

(<+l)m  /  1771  \  /  771 

22  (a(*+  1  ))/X/  =  5Z(a(l'))j^  22  akt 


;=.+i 


\j=* 

im  m 


fc=  1 


Setting  l  =  j  +  k  and  p  —  j  yields  i  <  p  <  im  and  1  <  /  —  p  =  fc  <  m,  and 

(t+l)m  (i+l)m  / /  — 1  \ 

52  (a(i  4-  i))<x'  =  51  52(a(l))Pa/-PM; 

(=«+i  f=«+i  \p= I  / 

Thus,  equating  coefficients  yields  the  recurrence 

/-i 

(a(i  4-  1)),  =  52(a(z))pa,-p 


p-  * 


Recurrence  for  (6(j)/  is  analogous. 
From  (2.1),  we  therefore  obtain 


m(n— 1) 


a?  =  coefficient  of  tk  in  (w{t)  52  a]  tfJ  +  22  eb(P(0)’U(0)-') 


j=i 


I +;=n 


=  12  12  eoMOU&U)), 

J=l  i+j=n  p+q=k 

In  particular,  a[  =  e10aic  4-  eoi^fc- 

For  an  integral  parametric  curve  r(t),  a  straightforward  computation  shows 
that  the  specialize  to 

1  <  it  <  n  -  I 


Qfc  =  < 


1 

°>k 
n—  t 


+  E.+;=n  'Lp+q=k  ei]{a(i))p{b(j))q  n  <  k  <  (n  -  1) 


m 


Ip+,=*  cf-,(a(i))p(6(j)), 


(n  —  l)m  <  i  <  nm 
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2.1.2  Derivation  of  the  Method 
2. 1.2.1  Rank  of  the  Linear  System 

Having  explained  how  to  obtain  the  a£,  we  now  show  that  the  coefficient  matrix 
of  the  linear  system  defined  by  setting  a£  =  0,  k  =  1,2, ... ,  nm,  has  rank  at  least 
<^(n).  We  are  able  to  determine  a  nontrivial  solution  to  unknown  coefficients  by 
setting  one  of  the  coefficient  to  1  and  solving  the  system  a"  =  0.  a2  =  0, . . . ,  a"  =  0 
for  s  >  <p(n)  chosen  such  that  the  rank  is  y >(n ). 

Let  en  =  (eio,  «oii  e2oi  eo2i  •  •  • .  eno,  e(„_!)  i, . . . ,  e!  (n_i),  e0n)r  be  the  vector 


of  unknowns,  and  write  the  system  of  equations  a"  =  0,  —  0 . a"m  =  0  in 

matrix  form: 

Amnen  =  0  (2.2) 


Note  that  Amn  is  a  nm  by  <j(n)  +  1  matrix.  Furthermore,  the  maximum  rank  of 
Amn  is  9(n)  +  1  since  m  >  n  and  nm  >  «p(n)  +  1.  Example  3.2  shows  matrix  A32 
symbolically. 

Example  2.2  For  m  =  3  and  n  =  2.  A32  is 


«iCo 

bic0 

0 

0 

0 

aiCi  +  a2co 

^lcl  +  62C<) 

O 

a^bi 

bj 

ajCo  +  a2Ci  +  a3Co 

b\C2  +  62C[  +  63cq 

2a  [a2 

Q{bn  +  026i 

26[62 

aic3  +  a2c2  +  a3Ci 

b\C3  +  62c2  +  63c  1 

2aia3  +  a* 

at63  +  a262  +  a36i 

26163  +  65 

a2c3  +  a3c2 

62c3  +  63c2 

2a2a3 

a263  +  a362 

26063 

a3c3 

b  3C3 

O 

a3 

a363 

bl 

When  computing  the  local  implicit  approximation  gn[x,y)  of  r(f),  if  the  rank 
of  Amn  is  at  least  ip(n)  then  we  can  select  one  coefficient  of  gn[x,y )  to  be  1  and 
determine  the  others  by  selecting  the  first  3  rows  of  (2.2)  and  choosing  s  such  that 
the  system  has  rank  >?(n) 
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Let  fm( x,  y)  =  0  be  the  exact  irreducible  implicit  form  of  r(<)  with  Fm(i,  y,  z) 
as  its  corresponding  homogeneous  form,  and  let  fn,  n  <  m,  be  the  degree  n  ini¬ 
tial  segment  of  fm  with  its  corresponding  homogeneous  form  Fn(x,y,x),  that  is, 
/m(x,y)  =  fn(x,y)+terms  with  degree  >  n.  Also,  let  b, V  =  Fn(p(t),q(t),w{t)) 
and  brrtn  —  (^li  ■  i  &mn) 

Lemma  2.1 

1.  /n(x,y)  is  the  zero  polynomial  if  and  only  if  bmn  =  0. 

2.  If  fn(x,y)  is  a  nonzero  polynomial  then  bi  =  62  =  •  •  •  =  6„  =  0  and  bmn  ^  0. 


27 


for  every  t.  By  comparison,  we  have  bx  =  frj  =  •  •  •  =  bn  =  0.  The  rest  of  part  2 
follows  from  part  1.  _ 

Since  /n(x,  y)  is  an  initial  segment  of  /m(x,  y),  it  could  be  either  a  zero  polyno¬ 
mial  or  a  nonzero  polynomial  with  zero  or  nonzero  bmn  vectors  respectively.  If  bmn 
is  known  beforehand,  the  coefficient  vector  en  of  fn(x,y)  is  uniquely  determined 
by  Amne„  =  bmn  which  is  an  overdetermined  linear  system.  Note  that,  for  a  fixed 
n,  the  elements  of  the  matrix  Amn  depend  only  on  the  coefficients  of  p(f),  q(t ),  and 
w(t).  The  following  results  characterize  the  rank  of  Amn. 

Lemma  2.2  If  r(t)  is  a  properly  parameterized  rational  plane  curve  of  degree  m 
then  for  n  <  m,  we  have 

rank(  Amn)  =  ^(n)  +  1 

Proof:  Suppose,  knowing  bmn,  we  want  to  determine  the  coefficients  of  fn{x.y)  by 
solving  the  overdetermined  linear  system  Amne„  =  bmn,  where  e„  is  the  coefficient 
vector  of  the  general  degree  n  polynomial.  Since  Amn  depends  on  the  coefficients 
of  r(f),  we  consider  the  following  two  cases. 

First,  if  r(t)  is  a  curve  such  that  the  corresponding  /n(x,  y)  is  the  zero  poly¬ 
nomial.  then  b„,n  =  0  and  Amnen  =  0  is  a  homogeneous  system.  If  rank(  Amn)  < 
p(n)  +  1.  there  will  be  infinitely  many  nontrivial  solutions  as  well  as  the  triv¬ 
ial  solution  for  this  linear  system.  This  cannot  be  true  by  Lemma  2.1.  Thus 
ranfc(Am„)  =  <p(n)  +  1. 

Second,  if  r(t)  is  a  curve  such  that  the  corresponding  /n(x,y)  is  a  nonzero 
polynomial,  then  bmn  ^  0  and  Amne„  =  bmn  is  a  consistent  non- homogeneous 
system  since  there  is  always  a  solution,  that  is  with  etJ  the  coefficients  of  fn{x,y). 
Suppose  rank(Amn)  <  9(n)+l.  this  system  will  have  infinitely  many  solutions.  Let 
e‘  be  one  of  the  infinitely  many  solutions  and  e"  #  e„,  where  e„  is  the  coefficient 
vector  of  /n(x.y).  Let  also  hn(x.y\  be  the  corresponding  polynomial  of  e*  and 

h^lx.y)  =  hn(x.y )  -  rms  <>/  j^'x.y)  with  degree  >  n 
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Let  Hm(r.y,z)  and  Hn{x,y,z)  be  the  homogeneous  polynomials  of  hrn(x,y)  and 
hn(x,y)  respectively.  Since  Amnen  =  Amne‘  = 

Hn{p{i),q(t),w{t))  =  Fn(p(t),  q(t),  w(t)) 


for  every  t,  and  thus 

Hm(p(t),q(t),w(t))  =  Fm(p{t),q(t),w{t))  =  0 

for  all  t.  We  then  have 

rmf  P(£)  ?(0  \  _  LmfP(0  <?(*)  v  _  q 

'  U}(t)'  w(t)  W'(t)'  Ut(f) 

for  all  but  finitely  many  f.  Hence  fm{x.y)  and  hm{x.y )  represent  the  same  alge¬ 
braic  curve.  Since  fn(x,y)  hn(x,y),  fm(x.y)  and  hm(x,y)  differ  by  more  than 
a  constant  factor  which  contradicts  that  the  equation  of  an  irreducible  curve  is 
unique  to  within  a  constant  factor.  Therefore  rank(Amn)  must  be  p(n)  +  1.  ■ 

Lemma  2.3  If  r(t)  is  a  properly  parameterized  rational  plane  curve  of  degree  m 
then  for  n  =  m.  we  have 

rank(  Amm)  =  p(m) 

Proof:  If  n  =  m,  fm(  p(t)/w(t),  q(t)/w(  t))  =0  for  all  t  except  finitely  many  t  where 
wit)  =  0.  and  then  Fm(p(t),  q{t),w\t))  =  0  for  all  t.  we  thus  have  Ammem  =  0 
which  is  an  overdetermined  linear  homogeneous  system.  Since  we  will  have  only 
trivial  solution  if  rank{ Amm)  =  'p(m)  +  1.  rank(Amm)  must  be  less  than  or  equal 
to  p{m I. 

Suppose  r  =  rankl Amm)  <  m ) ,  then  the  solution  space  of  the  overdeter¬ 
mined  homogeneous  system  has  as  basis  p  =  ^(m)  +  1  — r  linearly  independent 
vectors  and  every  solution  of  this  system  is  the  linear  combination  of  these  p  solu¬ 
tions.  Now  suppose  that  r  <  -p(m),  then  the  system  has  a  solution  space  spanned 
by  p  >  2  linearly  independent  vectors,  say  e^,  e;,,  •  •  ■  .  e*,.  Let  ff'ix.y)  be  the 
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corresponding  polynomial  with  coefficient  vector  e^  and  F™(x,y,z)  be  the  homo¬ 
geneous  form  of  fln{x,y  ),  i  =  1. 2, . . . ,  p.  Since 

FT(MM‘)M‘))  =  0  =  (i»«))"77(44-,  rz'l 

W{t)  w{t) 

for  all  t,  1  <  i  <  p,  we  have  f™(p(t)/w{t),q(t)/iv(t))  =  0  for  all  t  with  finitely  many 
exceptions  where  w(t)  =  0,  1  <  i  <  p.  Thus  the  irreducible  curve  /m(x,y)  =  0 
can  be  represented  by  f™{x,  y),i  =  1, 2, . . . ,  p,  which  are  not  different  within  only 
a  constant  factor  because  e^,  e^,  •  •  • ,  e^  are  linearly  independent.  By  the  above 
arguments,  we  can  conclude  that  rank(  Amm)  =  y’(m).  _ 

By  assigning  one  variable  to  be  1,  the  existence  of  a  nontrivial  solution  of 
Ammem  =  0  is  guaranteed  by  Lemma  2.3.  and  it  is  the  coefficient  vector  of  the 
exact  implicitization  of  r(t). 

Observe  that  Lemmas  2.2  and  2.3  are  not  valid  for  improperly  parameterized 
rational  plane  curves,  as  shown  in  the  following  example: 

Example  2.3  Let  C|(0  =  {x(t),y(t))  =  {i2  4-  2 1,  t 4  Hr  4 f3  4-  6t2  +  4 1).  Since  x(t)  = 
s(t)  and  y(t)  =  (s(t))2,  where  s(t)  =  t2  4-2,  c^f)  is  improperly  parameterized.  For 
n  =  2  the  rank  of  A42  is  4  whereas  ^(2)  4-  1  is  5.  g 

We  summarize  the  lemmas  above  in  the  following 

Theorem  2.1  If  r  (t)  is  a  properly  parameterized  rational  plane  curve  of  degree  m 


rank{  Amn)  = 


V?(  n )  4-  1  if  n  <  m 


[  9(n)  ifn  =  m 

where  n  is  the  degree  of  the  approximant,  ^p(n)  =  (n2  4-  3n  —  2)/2  is  the  number  of 
coefficients  on  which  the  approximant  depends. 


2. 1.2.2  The  Algorithm 

Because  of  Theorem  2.1,  we  may  compute  the  degree  n  local  implicit  approxi¬ 


mation  as  follows: 
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Let 

Ben  =  0  (2.3) 

be  the  subsystem  of  (2.2)  consisting  of  the  first  s  equations  of  Amnen  =  0  such 
that  B  has  rank  <p(n).  If  the  origin  is  not  a  singular  curve  point,  then  augment  the 
system  (2.3)  with  an  equation  e0i  =  1  or  e10  =  1  according  to  whether  x'(0)  ^  0 
or  y'(0)  7^  0.  If  the  origin  is  a  singular  curve  point,  on  the  other  hand,  then  (2.3) 
is  augmented  by  the  equation  e,j  =  1  where  the  indices  i  and  j  are  selected  by 


inspecting  the  system.  In  this  way,  a  linear  system 

Cen  =  b  (2.4) 

is  obtained  that  has  a  nontrivial  solution  for  en.  Since  a™  =  0,  a?  =  0 . a*  =  0. 


the  curve  gn  so  defined  must  have  contact  of  degree  at  least  s  to  r(<)  at  the  origin. 

There  may  be  cases  in  which  the  system  (2.4)  is  inconsistent,  i.e.,  the  augmented 
matrix  [C.  bj  is  of  rank  <p(n)  +2  while  C  has  rank  ip(n)  +  1.  In  this  case,  the  linear 
system  can  be  modified  to  ensure  consistency.  For  instance,  when  computing  g 2 
of  c 2(t)  =  ( t ,  '9),  q|  should  be  removed  from,  and  Qi8  =  0  should  be  added  to  the 

system  (2.4),  resulting  in  e01  -  1  =  0,  a?  =  0,  a\  =  0 . a|  =  0.  a?0  =  0.  a?8  =  0. 

In  this  way  a  g2(x.y)  =  xj  is  obtained,  an  approximation  that  has  eighth  order  of 
contact  and  is  irreducible. 

2. 1.2.3  Irreducibility  of  Implicit  Approximations 

When  the  origin  is  a  regular  curve  point  we  show  that  the  implicit  approxima¬ 
tion  gn(x,y)  of  r(t)  at  the  origin  is  irreducible  whenever  the  linear  system  (2.4)  is 
consistent.  Note  that  the  local  implicit  approximations  of  different  degrees  derived 
have  the  same  linear  terms  if  the  equations  augmented  to  (2.3)  are  the  same.  In 
the  following  lemmas,  we  assume  that  (2.  li  ;s  a  consistent  system.  Also,  let  s(n ) 
be  the  order  of  contact  made  by  the  degree  •>  unpiicit  approximation  gn(x.y). 


31 


Lemma  2.4  If  gn(x,y )  and  gn~l(x,y)  of  r(t)  at  the  nonsingular  point  (0,0)  are 
computed  by  augmenting  the  system  (2.3)  the  same  a  =  0,  where  a  =  0  is  either 
e10  —  1  =  0  or  e0i  —  1  =  0,  then  we  have  s(n  —  1)  <  s(n). 

Proof:  Let  gn{x,y)  =  %(gn{x,y)  +  gn~1{x,y)).  Suppose s(n- 1)  >  s(n),  the gn{x,y) 
so  defined  has  the  following  four  properties:  (1)  gn  is  of  degree  n.  (2)  gn  gn.  (3) 
gn  has  the  same  linear  terms  as  gn  since  gn  and  gn~l  have  the  identical  linear  terms. 
(4)  gn  has  the  order  of  contact  larger  than  or  equal  to  s(n)  since  s(n  —  1)  >  s(n) 
is  assumed.  From  properties  1,  3,  and  4,  the  coefficients  of  gn  satisfy  the  linear 
system  that  is  used  to  compute  the  coefficients  of  gn;  but  property  2  contradicts 
the  uniqueness  of  solution  of  a  nonsingular  linear  system.  Thus  s(n  —  1)  <  s(n).  ■ 
By  induction  and  Lemma  2.4,  we  can  show  that  s(n)  is  strictly  monotone. 

Lemma  2.5  At  the  nonsingular  point  (0,0),  the  degree  n  local  implicit  approx¬ 
imation  gn{x,y )  of  the  degree  m  >  n  properly  parameterized  parametric  curve 
r (t)  =  [x(t),y(t))  is  irreducible. 

Proof:  Suppose  gn(x,y)  is  reducible,  and  gn  =  gkgl,  where  n  =  jfc  +  /  and  k.l  > 
0.  Since  gn  contains  linear  terms,  one  of  the  gk  and  gl  must  have  a  constant 
term.  Let  gk{x,y)  =  Pi}x'yJ ,  where  pio  and  p0i  are  not  both  zero,  and 

9l(x,y)  =  EU;=o <U3x'y},  where  q00  #  0.  Let  also  that  gn(x(t),  y{t))  =  £^0,4', 

gk(x(t).  y{t))  =  Y.?=\Ptt\  and  p*(x(t),  y(t))  =  E^'0  7<t\  where  =  ?oo-  We  thus 

have 

mn  mk  ml 

£Q.*‘  =  (I>n(i:7,o 

t=l  1=1  i=0 

The  coefficients  of  gn{x,y)  are  computed  by  solving  the  nonsingular  system  as  (2.4) 
for  some  s  >  t p{n).  Moreover.  s(n),  the  order  of  contact  of  gn ,  is  greater  than  or 
equal  to  s.  Thus  the  coefficients  of  gn  satisfy  the  linear  system  a  =  0,q,  =  0,  aj  = 
0 . =  0.  where,  without  loss  of  generality,  we  assume  that  a  =  el0—  1  =  0. 
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The  above  linear  system  can  be  represented  in  terms  of  3,  and  7,  as  follows: 

qooPio  =  1 
4oo^i  =  0 
<?oo/?2  +  di7i  =  0 

< 

<700/^3  +  d27i  +  A  72  =  0 


1  ?00/3j(n)  +  A(n)-l7l  + - b  &7»(n)-l  =  0 

which  imphes  qooPio  =  1  and  3\  =  32  =  ■  ■  •  =  3,(n)  =  0.  Thus  gk  has  either  order 
of  contact  larger  than  or  equal  to  s(n)  if  s(n)  <  km,  or  gk(x(t),  y(t))  =  0  for  all  t 
if  s(n)  >  km.  The  first  result  contradicts  the  fact  that  s(n)  is  strictly  monotone, 
and  the  second  contradicts  the  irreducibility  of  the  exact  implicitization  of  r(f). 
Thus  gn  is  irreducible.  _ 

When  the  origin  is  a  singular  curve  point,  the  implicit  approximation  is  not 
always  irreducible.  For  example,  the  degree  n  implicit  approximation  of  c 3(t)  = 
(t2,<5),  with  implicit  form  i5  -  y2  =  0.  is  yn  =  0.  Note  that  y  =  0  is  the  curve 
tangent  at  the  origin. 

2.1.3  Error  Analysis 

2. 1.3.1  Quality  of  the  Approximation 

Given  e,  let  T(e,  n)  >  0  be  such  that  for  all  |<|  <  T(e,  n)  the  orthogonal  distance 
d{t,n)  between  point  (x(t),y(t))  and  thedegreen  approximation  gn{x.  y)  =  0  is  less 
than  e,  assuming  that  ( x(t),y{t ))  is  a  regular  curve  point  for  all  |i|  <  T{e,n).  The 
distance  d(tp,n)  from  a  point  P  =  (xp,yp)  =  [x{tp),  y(tp))  on  the  curve  r(t)  to  the 
degree  n  approximation  gn(x,y)  =  0  is  the  solution  of  a  difficult  nonlinear  system. 
A  reasonable  estimate  of  d(tp,n)  would  be  the  distance  to  the  gn(x,y)  =  0  in  a 
direction  orthogonal  to  the  level  curve  gn{x.y)  =  c,  where  c  =  gn{xp,yp),  denoted 
bv  d  ( tp ,  n).  Note  that  d'(t,  n)  >  d(t,n)  since  d(t,n)  is  the  shortest  distance  from 
the  point  to  the  curve.  Let  P'  =  {x'p,yp)  be  the  point  on  gn{x.y)  =  0  on  which 
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gn(x,y)  =  0  intersects  the  line  orthogonal  to  level  curve  gn{x.y)  =  c  at  P.  see 
Figure  2.1.  The  value  of  gn  at  P'  expressed  as  its  Taylor  series  about  P  is 


9n(xr  y'p )  =  9n(xp,  yP)  +  d'{tp,n)  ■  Vgn(xpj  yp)  +  higher  order  terms 


Taking  the  linear  term,  since  gn{ xp,y'p)  =  0.  d{tp,n)  can  be  approximated  by 
d"(tp,n)  where 


_  9  {xpiyP) _  9n(x(tp)i  y(tp)) 

|(Vgn(xp, yp)||  "  [{^(^(^),y(g))2  +  (<7yn(x(^),y(^)))2]1/2 


Note  that  d”(t,n)  may  be  less  than,  greater  than,  or  equal  to  d{t,n)  although 
d‘(t,n)  is  always  greater  or  equal  to  d(t,n). 

We  have  found  no  method  for  computing  T(e,n)  analytically.  However,  in 
practice  we  only  need  a  method  of  obtaining  a  reasonably  good  estimate  of  T(e,  n). 
Thus  it  is  desirable  to  determine  T'(e,n ),  for  given  e  and  n,  such  that  d“(t,n)  <  e 
for  |f|  <  T'(e,n). 

Since  2ab  <  a2  +  b2  for  any  positive  a  and  6,  we  have  \J\a\  jij  <  < 

so  that 


d  (<, n)  <  d(t.n)  = 


V2gn(x(t),y(t)) 


When  tracing  r(/),  we  can  detect  the  first  value  of  t  such  that  d(t,n)  <  t  and 
d{t  +  Al.n)  >  e.  where  is  the  step  distance  tor  t. 


2. 1.3.2  Curve  Translation  to  the  Origin 

In  the  derivation  of  the  approximant  we  assume  that  r(0)  =  (0,0),  i.e.  we 
require  that  r(<)  be  translated  to  the  origin  and  reparameterized.  Since  this  may 
incur  additional  inaccuracies  we  comment  on  it  now. 

Translation  of  r (t)  to  the  origin  is  a  simple  operation  that  incurs  only  small 
errors.  For,  with  p  =  (u,v)  as  the  curve  point  to  be  brought  to  the  origin,  the 
translated  curve  is  simply 

xi(t)  =  x(t)  —  u  =  I  p(t)  —  uw(t))/w[t) 

'Ji(t)  =  ’Jit)  ~  r  =  <7(0  -  eu:(0)/«tf(0 
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So,  we  have  to  subtract  two  polynomials  in  order  to  bring  p  to  the  origin. 

Now  assume  that  r(t0)  =  p,  and  consider  reparameterization  such  that  p  not 
only  is  moved  to  the  origin  but  that  also  to  =  0,  for  the  reparameterized  curve. 
Here  we  need  to  substitute  t  +  t0  for  t ,  i.e. 

x2(t)  =  Xi(t  +  t0) 

U2 (t)  =  !/i(t  4-  to) 

is  the  final  curve.  As  observed  in  the  introduction,  although  substitution  is  con¬ 
ceptually  simple,  it  nonetheless  may  introduce  numerical  errors  that  could  be  sig¬ 
nificant.  According  to  experiments  by  Prakash  and  Patrikalakis  [54],  Kahan's 
method  described  in  [39]  exhibits  good  numerical  stability,  and  offers  one  method 
of  implementing  the  needed  reparameterization. 

A  second  method  would  be  to  avoid  reparameterization  altogether  by  reformu¬ 
lating  the  derivation  of  the  approximant  given  before.  That  is.  we  consider  rif) 
containing  the  origin  at  which  t  is  not  necessarily  0.  seeking  again  an  implicit  ap¬ 
proximant  at  the  origin.  Clearly  this  is  possible  and  requires  only  straightforward 
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modifications  of  our  method.  In  fact,  even  translation  of  the  point  to  the  origin 
can  be  avoided  by  such  modifications.  The  details  are  routine. 

2.1.4  Experiments 

A  good  local  approximation  of  a  curve  will  provide  a  more  accurate  local 
approximation  and  larger  interval  T(e,n)  of  approximation,  for  a  given  e,  when 
the  degree  n  of  the  approximation  increases.  The  local  implicit  approximation 
gn(x,y)  =  0  of  r (t)  is  determined  by  y?(n)  linear  conditions  imposed  on  its  coeffi¬ 
cients,  where  <j?(n)  is  the  degree  of  freedom  of  gn(x,y).  Thus,  as  n  increases,  more 
conditions  can  be  satisfied  and  finally  the  exact  impiicitization  is  obtained  when 
n  =  m.  Hence  our  local  implicit  approximation  is  capable  of  approximating  a  given 
curve  not  only  locally  but  also  globally  in  the  sense  that  T(e,n),  for  a  given  c,  will 
be  larger  when  n  increases.  On  the  other  hand,  a  local  explicit  approximation  is 
limited  due  to  the  asymmetry  introduced  by  making  one  variable  an  explicit  func¬ 
tion  of  the  other.  Thus,  an  local  explicit  approximation  can  only  approximates  the 
given  curve  locally  for  |x|  <  R,  where  R  is  its  radius  of  convergence,  no  matter 
how  high  the  degree  of  approximant  is. 

We  give  as  example  the  approximation  of  several  parametric  curves  that  are 
not  singular  at  the  origin,  showing  both  implicit  and  explicit  approximations. 

Example  2.4  Four  curve  examples  are  shown  below. 

c4(  t)  =  ( t6  +  t'°  -  2 f3  +  3 12  +  12 1. 16  -  t'°  +  t4  -  4f3  -  2 t2  +  24<) 

c s(t)  =  (3 16  -  4<5  -  8 f3  +  6 t2  4-  3 1,  -3 16  +  4 15  +  ot*  -  6f3  -  St2  +  3 1) 

Cg(  t)  =  (3 16  +  ts-  2 14  +  3Sf3  -  ot2  -  14 f,  t6  -  12 15  -  2 14  +  2 f3  -  7f2  +  13f) 
cT(0  =  ((t6  +  3t5  —6t4  +  4t3  —  36t2  +36t)/w(t),  (3t6  +  t5  —  2t‘i  +  39t3  —  69t2  +  33t)/w( t)) 
where  w(t)  =  7 16  +  10t5  -I-  9<4  -I-  6t2  +  3t  +  7. 

The  curves  of  c4(f).  c 5(t),  c $(t)  and  Cr(t)  with  t  in  [—1,1],  and  their  lo¬ 
cal  implicit  approximations  and  local  explicit  approximation  are  shown  in  Fig¬ 
ures  2.2.  2.3.  2.4  and  2.5.  Note  the  good  quality  of  local  implicit  approximation. 
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Tables  2.1  2.2  and  2.3,  for  c4(f), c5(f)  and  c «(£)  respectively,  list  the  y-values  of  a 
sequence  of  x-values  to  quantify  how  accurately  the  low  degree  local  implicit  forms 
approximate  the  original  curves.  The  corresponding  values  of  local  explicit  forms 
are  also  listed  for  comparison.  For  such  examples,  we  observe  that 

(a)  For  local  implicit  approximation.  d(£,  n  +  k)  <  d(£,  n)  for  t  in  [—  1. 1]  and  k  >  1. 
In  addition,  T(e,n)  <  T{e,n  +  k)  for  k  >  1. 

(b)  r(e,  2)  of  local  implicit  approximation  is  greater  than  T(e,6)  of  local  explicit 
approximation. 

(c)  Degree  2  and  3  local  implicit  approximations  give  very  accurate  approximations 
on  a  reasonable  range  of  t. 

(d)  Degree  5  local  implicit  approximation  approximates  the  original  curve  very 

precisely  at  least  for  —  1  <  t  <  1.  ■ 

When  computing  an  explicit  approximation  y  =  h(x)  directly  from  the  curve  r|f), 
we  first  compute  the  degree  n  power  series  t  =  £*=1  dxx'  from  x  =  x(t),  and  then 
substitute  it  for  t  in  y  =  y(t).  As  a  result,  only  the  first  n  coefficients  of  h{x) 
are  exact  and  the  remaining  coefficients  obtained  in  the  computation  should  be 
discarded.  Moreover,  substituting  t  =  HILj  f°r  f  in  y  =  y( t)  is  not  a  cheap 
computation,  especially  for  high  degree  local  explicit  approximations.  Hence,  the 
computation  of  a  local  explicit  approximations  of  a  parametric  curve  directly  from 
r|  t)  is  more  costly  than  the  implicit  form.  In  general,  the  computation  of  local 
implicit  approximation  involves  generating  the  a,"*  and  solving  the  linear  system, 
which  is  fairly  efficient  for  low-degree  approximation. 

The  local  explicit  approximation  is  an  analytic  function  that  does  not  exist  at 
a  curve  singularity.  In  contrast,  a  local  implicit  approximation  always  exists. 

Example  2.5  Local  implicit  approximations  can  be  derived  at  singularities,  in¬ 
cluding  cusps,  where  local  explicit  approximation  fails.  Let  cs(<)  =  (of3  +  2 £*.  t*  — 
5f3  +  2f'2)  with  the  implicit  form /4(j. y j  =  — .r 1  r5oxJ -r6S3x*y  +  1325xy* 4-625yJ  — 
336i*  -r  672xy  -  336 y*.  The  origin  is  a  cusp  of  cs(£)  with  tangent  x  -  y  =  0.  The 


aegree  1  locai  implicit  approximation- 


Figure  2.5  Curve  c-(t) 
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degree  2  local  implicit  approximation  is  a  double  line  (x  —  y )2,  which  is  the  best 
degree  2  approximation  one  can  derive  at  cusp.  The  degree  3  local  implicit  approx¬ 
imation  is  x2  -  2 xy+y2  -  0. 16259766x3  -  2.0356445x2y  -  3.940918xy2  - 1 .8608398y3 
which  shows  very  nice  approximation  to  the  c8(2)  with  2  in  [—1, 1],  see  Figure  2.6. 
As  a  next  example,  we  consider  c9(2)  =  ((52s  —  162‘*-f  1023-|-422)/u;(2),  (ts  +  t4  +  2t3 — 
16t2)/w(t)),  where  w(t)  =  0.123  4-  0.122  —  22  +  12.5.  The  c9(2)  is  a  singular  curve 
with  a  cusp  at  the  origin  and  a  self-intersection  as  well,  as  shown  in  Figure  2.7. 
Figure  2.8  shows  the  degree  3  and  degree  4  local  implicit  approximations  of  c9(2). 
The  degree  4  local  implicit  approximation  shows  remarkable  performance.  ■ 


2.2  Local  Implicit  Approximation  of  Parametric  Surfaces 


We  derive  an  implicit  surface  g(x,y,z)  =  0  that  approximates  the  parametric 
surface  P(s,  2)  =  (x(s,  2),  y(s,  2),  z(s ,  2))  at  the  origin  to  a  specified  order  of  contact, 
using  the  method  of  Section  2.1.  Let 


P(s,2)  =  (x(s,2),y(s,  2),z(s,2))  =  ( 


Pi3’  0  <?(M)  r(s, 2) 

w(s,t )’  w(s,  2)  ’  w(s,t) 


be  a  rational  parametric  surface  of  total  degree  m  containing  the  origin,  where 


p(M)  =  XI  a>js'tJ-  q{s,t)=  £3  bt]s'tJ, 

>+j=i  i+j=i 


m 

r(s,t)=  53  cv3'tJ’  w(s,t)=  J3  dijs't2 

•+J  =  1  t+;=0 

with  a,j,  bi},  dj,  i+j  =  m,  not  all  zeros  and  doo  ^  0.  It  has  been  shown  by  Macaulay 
[46]  that  a  parametric  surface  of  the  above  form  has  an  irreducible  implicit  iorm 
/a(x.  y,  -)  =  0  of  degree  d  <  m2. 

Let  gn{x,  y,  z)  =  'E?+J+k=\  eijkXly} zk .  n  <  m2.  with  symbolic  coefficients  el]k 
be  a  general  implicit  form  of  degree  n  surface  that  contains  the  origin.  Since 
gn(x.y,z )  =  0  is  unique  up  to  a  constant  factor,  it  has  degrees  of  freedom  p(n)  = 
(fn  +  l)(n  +  2)(n  +  3))/6-2. 


4 


(t)  :  degree  5  rational  curve, 

he  origin  is  a  cusp  and  a  self— intersection  point, 
abel  i  :  cegree  i  local  implicit  approximation. 


4-2  y 
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When  substituting  x(s,  t),  y(s,  t)  and  z(s,  t)  into  gn(x ,  y,  2),  we  obtain 

n/_z  ^  ,  iXX  Gn(p(s,t),q(s,t),r{S't),w(s,t))  £™=i 

g  {x{s,t),y{s,t),z(s,t))  =  -  - - -  —l——- 

(w(s.t))n  (w(s,t))n 

where  Gn(x,y,  z,w)  is  the  homogeneous  form  of  gn(x,y,z),  and  the  a,j  are  linear 
combinations  of  e.jfc.  The  local  implicit  approximation  gn(x,  y ,  2)  of  the  parametric 
surface  P(s,  t)  is  computed  as  in  the  case  of  the  curve  approximation.  The  following 
section  shows  a  recursive  derivation  of  a,,  that  obviates  the  need  to  explicitly 
substitute. 

2.2.1  A  Recurrence  for  atJ 

Let  Gn{x,y,z,w)  and  Gn~l(x,y,  z,w)  denote  the  homogeneous  polynomials  of 
gn(x,y,z)  and  gn~1(x,y,  2),  respectively.  Since 


we  have 


gn(x,y,z)  =  gn~l(x,y,z)  +  £  e.^x't/V 

>+j+k=n 


Gn(x,y,z,w)  =  wGn~l{x,y,z,w)  +  ^  eijkx'yJzk 

i+j+k=n 


We  define  (a(k))ij,  ( b(k))ij  and  ( c( A: ) ) , ^  by  setting 


and  similarly  for 


(p(s,t))k  =  (  £  avs'tJ)k=  £  («(*)), •>*’* 

<+.7=1  i+j=k 


{q{s,t))k  =  J2  (&(£)), 

1 +]=k 


Since 


»+;=* 

Arm  /  (/s—  Dm  \ 

(a(fc))I_,s'fJ  =  I  £  (a(fc  -  l))^^  ]p(s,f) 

l+j  =  Jc  \p+q=k- 1  / 
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as  shown  in  [48],  the  (a(A:))u  is  recursively  defined  as 


1 

0 

aH 


if  k  =  i  =  j  =  0 
if  k  =  0  and  i  +  j  >  1 
if  k  =  1  and  i  4-  j  >  1 
if  k  >  2  and  z  +  j  >  A: 


(*(*))<i  and  (c(*)).j  are  defined  analogously. 

Recursive  derivations  for  (a(fc)),-j,  {b(k))t]  and  (c(Ar))IJ  can  be  found  in  [4S]. 
Let  a?  and  o'1"1  be  the  coefficient  of  s't 3  in  G'ipis,  t),  q(s,t),r(s,  t),  w(s,  t))  and 
Gn_1(p(s,t),9(s,t),r(s,f),u;(s,t)),  respectively.  From  (2.5),  a”  can  be  derived 
from  the  a^f1 ,  where  1  <  k  <  i  and  1  <  /  <  j,  as  shown  in  the  following  formula: 


m(n— 1) 

=  coefficient  of  s' V  in  (w(s,  t)  ^  ct"~ls'tJ 

1+7=1 

+  L  «W,0K«,  0)*‘  (?(*,  t))k'(r(s ,  f))*3) 

fct+fcj+*3=n 
i  J 

=  E  E 

/,  =  l  h  =  l 

d"  5Z  5Z  e fcl  *2*3  (  a(  ^-1  )  )pi  <71  (  ^(  ^-2  )  )p^<72  (C(  ^3)  ) P3<J3 

*l+^2+*3=n  Pl+P2  +  P3  =  * 

71+73+93—7 

Note  that  a]}  —  +  eoio^ij  4-  eooici;. 

For  an  integral  parametric  surface  P(s.t)>  since  (a( Ar))tJ  =  0 .(b(k))tJ  =  0  and 
(c(fc))t>  =  0  for  1  +  j  <  k  and  a £  =  0  for  i  4-  j  >  (n  -  l)m.  we  have 


for  1  <  i  4-  j  <  n  —  1, 


for  n  <  i  +  j  <  (n  —  1  )m. 


Q1 7  aiJ  d"  ZZ  ZZ  e^l  *2*3  (  G(  ^'1  ))pi  71  (M  ^2  )  )p392  (C(  ^'3)  )p393 

*1  +*2  +  fc3  =  TX  Pl+Pl  +  P3=l 

71  +72  +  73  =  7 
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for  (n  —  l)m  <  i  +  j  <  nm , 


°L\)  ^  I  y  -■  (  a(  )  )pi9l  (^(^'2  )  )pj?2  (*-(  ^-3))p3l}3 

*:i+fcj+*3=n  PI  +  P2+P3=> 

<?l+«+<73=J 


2.2.2  Derivation  of  the  Method 


2. 2. 2.1  Rank  of  the  Linear  System 


Having  derived  q|*  ,  i+;  =  1,2,...,  nm.  for  the  degree  n  implicit  approximation 
gn(x,y,z),  we  write  the  system  of  linear  equations  a"0  =  0,  =  0,  c^o  = 

(nm)0  ~  a(nm— 1)1  —  ^ . ai(nm-l)  =  <^0(nm)  = 


0,  aj*!  =  0.  Qq2  =  0. 
0  in  matrix  form 


a 


Arnn^n  —  0 


(2.6) 


where  en  —  (eioo^oio,  ^ooi,  ^200^  ^lioi  eioi.  e02o*  eon,  ^002, . . . ,  eoon)1*  is  the  vector  of 
unknowns.  Amn  so  defined  is  ofdimension  ((nm  +  l)(nm+2))/2-l  by  p(n)  +  l,  and 
has  rank  at  most  p(n)  +  1.  As  in  the  curve  case,  the  rank  of  Amn  is  critical  when 
solving  for  the  unknown  coefficients  etJk-  The  following  theorem  characterizes  the 
rank  of  Amn. 


Theorem  2.2  IfP(s,  t)  is  a  properly  parameterized,  rational  surface  of  total  degree 
m.  then 


rank(Amn) 


p(n)  +  1  if  n  <  d 
p(n)  if  n  =  d 


where  d  (<  m2)  is  the  degree  of  the  implicit  form  ofP(s.t),  n  is  the  degree  of  the 
approzimant,  and  p(n)  is  the  degree  of  freedom  of  the  approximant. 


Proof:  Similar  to  proofs  of  Lemma  2.2  and  Lemma  2.3.  g 

As  a  result  of  Theorem  2.2,  it  is  clear  that  the  exact  implicitization  of  P(s.t) 
is  the  solution  of  Ammem  =  0.  with  one  variable  -?et  to  a  fixed  value. 
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2. 2. 2. 2  The  Algorithm 

We  compute  the  degree  n  implicit  approximation  gn{x,y,z)  as  follows:  Let 

Ben  =  0  (2.7) 

be  the  subsystem  of  (2.6)  that  consists  of  the  first  s  equation  of  (2.6)  such  that  B 
has  rank  p(n).  Augment  the  system  (2.7)  with  a  =  0,  where  a  is  determined  as 
follows:  if  the  origin  is  a  regular  surface  point,  a  is  eioo  —  1.  e0io  —  1.  or  eooi  —  1 
depending  on  the  gradient  of  the  surface  at  the  origin.  If  the  origin  is  a  singular 
surface  point,  then  a  =  et]k  —  1  where  the  indices  i,j,  and  k  are  selected  by 
inspection.  Thus,  a  linear  system 

Cen  =  b  (2.S) 

is  obtained.  System  (2.8)  may  be  an  inconsistent  system.  If  this  happens,  some 
equations  must  be  removed  from  (2.8)  to  ensure  the  consistency. 

One  alternative  for  handling  inconsistencies  is  that  we  replace  a  =  0  in  (2.8) 
with  en0o  —  1  =  0,  e0„o-  1  =  0,  or  eoon  -1=0  and  then  solve  it  as  usual.  Examples 
show  that  gn(x,y ,  z)  computed  by  this  method  can  be  of  the  form  (ax+by+cz)n  for 
some  a.  b.  c,  that  is,  it  degenerates  to  the  tangent  plane.  To  remove  this  degeneracy, 
we  do  the  following: 

1.  Solve  for  gl(x,y,z ),  and  compute 

2  ~  (91(*(s-0’y(M),-(s,*)))n 

t+7=n 

2.  Consider  the  linear  system  that  consists  of  the  first  s  equations  of  (2.6)  and 
a  =  0.  where  a  =  0  is  en00  —  1=0.  e0no  —  1  =  0,  or  eoon  —  1=0  and  s  is 
chosen  such  that  the  coefficient  matrix  of  the  system  has  rank  p(n). 

3.  Find  a  J(j  which  is  nonzero  and  augment  the  corresponding  a,;  =  0  to  the 
above  svstem.  then  solve  it. 
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This  computation  of  the  local  implicit  approximation  results  in  an  approximant 
that  has  roughly  n3/2-th  order  of  contact.  Thus,  when  raising  the  degree  of  the 
approximant,  the  order  of  contact  with  the  surface  P(s,  t)  grows  subquadratically. 

Example  2.6  Consider 

P(s,t)  =  {x{s,t)/w{s,t),y(s,t)/w{s,t),z(s,  t)/w(s,t)) 

where 


x(s,f)  =  — 200f2  +  12st  +  400*  -  200s2  -  10s 
y{s,t)  =  15<2  —  14st  +  lOt  —  11s2  +  400s 
z(s,t)  =  200t2  +  list  —  t  +  200s2  +  2s 
u;(s,  t)  =  100f2  -  200f  +  100s2  +  200 

We  compute  degree  2  and  degree  3  local  implicit  approximations 
g2{T,y,z)  = 

-108.44294z2  -  10.26463Syz  -  13.162097xz  +  381.19047z 

-  95.092S36y2  -  5.241114xy  -  l.SS09524y  -  94.85476x2  +  x 

and 

y3(x,y,z)  = 

1.3012126z3  +  5.16125yz2  -  46.6908lxz2  -  103.SlS164z2 

-  1.1158845y2z  +  15.62259Sxyz  4-  4.51S084yz 

+  1.267575x2z  +  lS0.40466xz  4-  381.19047*  -  3.6884814y3 

-  48.00386xy2  -  95.16589y2  +  5.56l3696x2y  -  6.1573525xy 

-  l.SS09524y  -  44.977395x3  -  94.34699x2  +  x 

Note  that  the  normal  of  /4(x.y,z).  the  exact  implicit  form  of  P(s.  t),  at  the 
origin  is  almost  parallel  to  z-axis.  Thus,  to  show  the  performance  of  the  local 
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Table  2.4  Maximal  deviations  between  the  intersection  curves 


li  n  :  local  implicit  form  of  degree  n.  le  n  :  local  explicit  form  of  degree  n. 


degree 

r=0.25 

r=0.50 

r=0.75 

r=1.25 

r=1.5Q 

r=1.75 

0.000000 

0.000004 

0.000035 

0.000837 

0.004964 

0.019612 

■ 

0.000289 

0.002531 

0.009630 

0.066627 

implicit  approximation,  we  intersect  the  cylinder  h(x,y,z)  =  x2  +  y2  —  r2  =  0 
with  the  surfaces  f4,g2,  and  g 3,  and  plot  the  intersection  curves  of  f4  =  0  n  h  = 
0,p3  =  0Pl/i  =  0,  and  g2  =  0  H  h  =  0  in  one  figure.  Figures  2.9  and  2.10,  show  the 
intersection  curves  in  cylindrical  coordinates,  for  r  =  0.25,0.5,0.75, 1.00, 1.25, 1.50, 
and  1.75.  Table  2.4  lists  the  maximal  deviations  between  the  intersection  curves 
f4  =  0  n  h  =  0  and  g3  =  0  D  h  =  0,  and  f4  =  0  n  h  =  0  and  j2=0fl/i  =  0.  ■ 


2.3  Remarks  on  Resultants 

Different  resultants  are  formulated  in  the  classical  literature  for  the  purpose 
of  eliminating  variables  from  systems  of  algebraic  equations.  Early  expositions  of 
several  formulations  are  found  in  Netto's  book  cited  in  the  references.  In  essence, 
resultants  constitute  a  projection.  A  well-known  problem  of  elimination  based  on 
resultants  is  the  possibility  of  obtaining  extraneous  factors.  For  example,  when 
implicitizing  the  parametric  sphere 

x(s,t)  =  (1  -  s2  _  t2)/{l  +  s’2  +  t2) 

y(s,t )  =  2t/[  1  +  s2  +  t2) 

z(s,  t)  —  2s/(l-f-s*-f-t*) 
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the  Sylvester  resultant  yields 

256(x  +  1)V  +y'  +  z2  -  l)2 
and  the  Dixon’s  resultant  yields 

— 64(x  +  l)(x2  +  y2  +  z2  -  1) 

In  each  case,  an  extraneous  factor  i  +  1,  to  some  power,  is  present. 

Technically,  a  resultant  is  based  on  formulating  a  system  of  linear  equations 
with  symbolic  coefficients.  This  is  especially  apparent  in  the  derivation  of  Sylvesters 
resultant.  Macaulay  recognized  that  extraneous  factors  technically  are  related  to 
dependent  equations,  and  that  they  can  be  eliminated  by  division  by  a  suitable 
minor  [47],  Modern  work  on  the  multivariate  resultant  tries  to  find  this  minor  algo¬ 
rithmically,  i.e.,  to  recognize  and  eliminate  extraneous  factors.  See,  e.g.,  [20,  13]. 
In  our  approach,  a  linear  system  is  formulated  numerically,  hence  dependencies 
among  the  equations  are  easy  to  recognize.  If  the  approximant  is  formulated  with 
the  exact  degree  of  the  implicit  form,  then  our  approach  determines  the  implicit 
form  without  extraneous  factors.  If  an  approximant  of  higher  degree  is  determined 
with  our  approach,  then  a  reducible  implicit  form  with  extraneous  factors  could 
be  generated. 
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3.  LOCAL  APPROXIMATIONS  OF  2-D  SURFACES 


In  solid  modeling  with  curved  surfaces,  a  number  of  desirable  surface  opera¬ 
tions,  including  offsetting,  spherical  blending  and  the  formation  of  Voronoi  sur¬ 
faces,  raise  difficult  mathematical  problems  that  must  be  solved  in  order  to  repre¬ 
sent  and  interrogate  the  resulting  surfaces.  Such  surfaces  cannot  be  easily  defined 
in  the  conventional  3-D  space.  In  contrast,  they  can  be  formulated  mathemati¬ 
cally,  with  the  theory  of  envelopes,  in  higher  dimensional  space  in  a  straightforward 
manner.  For  example,  given  the  surface  h(x,y,z)  =  0.  we  formulate  the  r-offset 
of  h  as  the  envelope  of  a  family  of  spheres  with  radius  r  whose  centers  lie  on  the 
surface  h  —  0.  This  formulation  results  in  a  system  of  four  polynomial  equations 
in  6-D  space.  Such  surfaces  are  generally  2-surfaces  in  Rn,  n  >  3,  and  are  defined 
by 

/i(xj,x2 . xn)  =  0 

/2(xi,x2 . xn)  =  0 

(3.1) 


f n—2  ( Xj ,  X2, - Xn)  =  0 

where  the  /,  are  polynomials  or  in  matrix  form  as  F(x)  =  0.  Although  the  exact 
closed-form  representation  of  such  a  surface 


/(x  i,x2,x3)  =  0  (3.2) 

could  be  derived  in  principle  by  elimination  methods  such  as  Grobner  basis  or 
resultant  techniques,  it  is  often  not  feasible  to  do  so  in  practice  due  to  the  high 
complexity  of  these  methods. 

For  surface  representations  i:i  h:gh-dsmensional  space  to  be  practical,  good 
algorithms  for  some  basic  operations  iuvc  to  be  developed.  Among  them,  as 
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mentioned  in  [32],  are  iocating  points  on  surface  intersections,  tracing  surface  in¬ 
tersections,  and  local  surface  approximations.  It  is  our  objective  here,  for  a  given 
surface  representation  F(x)  =  0  and  a  regular  point  x°  =  (x°,  x2, . . . ,  x°)  on  it.  to 
develop  computational  schemes  for  the  local  geometry  of  /(xi,x2,x3)  =  0  and  to 
investigate  approximation  schemes  which  derive  a  surface  of  lower  degree  that  ap¬ 
proximates  the  given  surface  locally  at  x°  =  (x°,  xi^x®)  in  (xt,  x2,  x3)-space.  These 
techniques  will  be  of  practical  interest  if  they  are  efficient  and  can  be  incorporated 
into  algorithms  for  surface  interrogation  such  as  computing  surface  intersections. 

In  this  chapter,  we  proceed  as  follows:  Section  1  presents  techniques  that  de¬ 
scribe  the  local  geometry  of  /(xi,x2,x3)  =  0.  Section  2  presents  an  algorithm  that 
computes  degree  two  implicit  approximants  of  /(x!,x2,x3)  =  0  without  actually 
computing  /.  Section  3  describes  the  computations  for  local  explicit  approximants 
to  /(xi,x2,x3)  =  0.  In  section  4,  we  derive  a  procedure  for  computing  parametric 
approximants. 

3.1  Local  Geometry  of  the  Projected  Surface 

For  a  given  2-surface  Sf  in  Rn  and  a  regular  surface  point  x°.  we  first  derive 
some  schemes  that  determine  the  normal  vector  and  tangent  vectors  of  S/  at  a 
regular  projected  point  x°  from  the  normals  and  tangents  of  Sf  at  x°.  We  then 
derive  an  algorithm  to  compute  the  normal  curvature  of  surface  5/  at  x°  in  a 
tangent  direction  v,  which  requires  only  the  information  provided  by  F(x)  =  0  at 
x°. 

In  the  following,  we  assume  that  the  surface  point  x°  and  its  projection  x°  are 
nonsingular.  Let  Txo(Sp)  denote  the  affine  tangent  space  to  Sp  at  x°.  that  is,  the 
set  of  all  tangent  vectors  to  surface  Sp  at  x°.  It  is  evident  that  T’xo(5/r)  is  the  null 
space  of  D Ffx°)  h  =  0.  and  that  Txo{Sp)  is  a  vector  space  of  dimension  2  since 
DF(x°l  has  rank  n  —  2.  Note  that  at  nonsinguiar  surface  points  the  tangent  space 


59 


of  the  surface  5/  at  x°,  denoted  as  7*o(S/),  has  the  same  dimension  as  Txo(Sp )  . 
That  is,  the  dimensionality  of  the  tangent  space  is  then  invariant  under  projection. 

A  vector  n  is  a  normal  to  surface  Sf  at  x°  if  n  •  h  =  0  for  every  h  G  Txo(Sf). 
The  normal  vectors  form  a  vector  space  of  dimension  n  —  2.  Since  D F(x°)  has 
maximal  rank  n  —  2,  the  gradient  vectors  V/i(x°), ....  V/n_2(x°)  form  a  basis  for 
the  normal  vector  space  of  surface  Sp  at  x°.  Thus  any  linear  combination  of  the 
gradients  vectors  is  a  vector  in  the  normal  space.  Note  that  the  normal  space  of 
surface  Sf  at  a  nonsingular  surface  point  is  of  dimension  1. 

3.1.1  Normal  Vector 

One  may  expect  that  computing  the  normal  vector  of  the  projected  surface  5/ 
at  point  x°  =  (xj,  x^x^)  from  the  D F(x°)  is  as  complex  as  the  elimination  process 
from  F(x)  =  0  to  /(xt,x2,x3)  =  0.  Indeed,  it  is  true  unless  we  approach  the 
computation  differently.  Instead  of  considering  the  global  surface  Sf ,  we  consider 
its  tangent  space  at  x°  and  the  tangent  plane  of  Sf  at  x°.  The  tangent  space  of 
Sp  at  x°  is  the  null  space  of 

£>F(x°)  (x  —  x°)  =  0  (3.3) 

from  which  the  equation  of  the  tangent  plane  of  5/  at  x° 

Ax\  +  Bi2  +  Cij  +  D  =  0  (3.4) 

can  be  obtained  by  linear  algebra  techniques  such  as  Gaussian  elimination  [43]. 
Considering  the  elimination  process,  we  obtain  the  following  algorithm  for  com¬ 
puting  the  normal  vector  of  5/  at  point  x°: 

Algorithm  3.1 

1.  Consider  C  =  Z7JV2  a, V/.fx0).  where  component  C,  of  C  is  a  the  linear 
combination  of  the  unknowns  o,.  i  =  1 . n  —  2. 
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2.  Solve  the  linear  system 

0+  =  0,  C%  =0,  . . . ,  Cn  =  0 

which  is  of  dimension  (n  —  3)  x  (n  —  2)  and  has  one  degree  of  freedom. 

3.  Let  (dj,  d2, . . . ,  d„_2)  be  a  solution  of  the  above  system  and  C  be  the  corre¬ 
sponding  C  in  step  1.  Then  the  first  three  components  of  C  are  the  normal 
vector  of  5/  at  the  regular  point  x°. 

Theorem  3.1  Algorithm  3.1  computes  the  normal  vector  of  the  projected  surface 
Sf  at  x°  =  (.11,22,13)  if  x°  is  a  nonsingular  point  on  S j . 

Proof:  Let  /  be  the  ideal  generated  by  the  system  (3.3).  Then 

C(x)  =  £a,V/,(x°)(x-x0) 

1=1 

lies  in  the  ideal  I.  Since  (3.3)  is  a  linear  system,  and  the  tangent  plane  (3.4)  of 
Sf  at  x°  can  be  obtained  by  eliminating  the  variables  x4,  x5, . . . ,  xn,  step  2  of  the 

algorithm  is  the  elimination  of  the  variables  24,x5, _ xn.  Therefore,  the  C(x) 

with  the  computed  coefficients  (di,  d2, . . . ,  d„_2)  is  in  /  n  A'fxx,  x2,  x3].  Hence  the 
vector  consisting  of  the  first  three  components  of  C  is  the  gradient  of  C(x)  which 
is  the  normal  vector  of  5/  at  x°  corresponding  to  the  normal  space  of  S p  at  x°.  ■ 
Note  that  this  computation  is  valid  only  when  x°  is  nonsingular.  Suppose  the 
computation  is  valid  at  the  singular  point  x°.  Then  the  matrix  after  applying 
Gaussian  elimination  to  £>F(x°)  will  not  have  maximal  rank,  which  contradicts  to 
the  nonsingularity  assumption  of  x°. 

3.1.2  Tangent  Vectors 

Since  the  affine  tangent  space  Txo(Sf\  <>t  at  a  nonsingular  surface  point  x° 
has  dimension  two  and  is  the  null  space  of  / 1 F '  x1 .  h  =  0.  we  let  ht  and  h2  be 
a  basis  of  Txo(^f)-  The  following  theorem  -:iuws  -:iat  tin*  corresponding  basis  of 
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Tx°(Sf)  are  the  natural  projections  of  hi  and  I12,  denoted  as  hi  and  h2  respectively, 
provided  that  hi  and  h2  are  linearly  independent. 

Theorem  3.2  If  hi  and  h2  nre  the  natural  projections  of  hi  and  h2,  respectively, 
into  (x\,xi,x$)-space  and  are  linearly  independent,  then  hi  and  Ji2  form  a  basis  of 

Tv(S,). 

Proof:  Since  DF(x°)  has  maximal  rank  n  —  2,  its  row-echelon  normal  form  is 


A 

B 

C 

0 

0 

0  • 

•  0 

X 

X 

X 

* 

0 

0  • 

•  0 

X 

X 

X 

X 

* 

0  • 

••  0 

X 

X 

X 

X 

X 

x  • 

*  •  * 

where  .4,  B ,  and  C  are  coefficients  of  linear  terms  in  (3.4),  *’s  are  nonzero  numbers 
and  the  x!s  represent  numbers  that  may  or  may  not  be  zero.  Since  ht  and  h2  form 
a  basis  of  Txo(Sf)  h\  and  h2  are  linearly  independent  and  satisfy 

(A,  B,  C,  0, 0, 0,  •  ■  • ,  0]  h  =  0 

Thus  ht  and  h2  satisfy 

[A.B.C]  h  =  0 

Furthermore,  hi  and  h2  are  linearly  independent.  Hence  hi  and  h2  form  a  basis 
of  7*o(67).  ■ 

For  a  given  tangent  h  to  Sf  at  a  nonsingular  surface  point  x°  the  tangent  to 
5/  at  x°  is  the  natural  projection  of  h  assuming  x°  is  nonsingular. 

Example  3.1  Consider  the  general  parametric  surface 

.r  =  M-M) 

y  =  h2{s,t ) 

Z  =  ki(s,t) 
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and  a  surface  point  (x°,y°,z0)  =  (hx(s0,  t°),  h2(s°,  i°),  h3(s°,  t0)).  In  five  dimen¬ 
sional  space,  we  write  the  surface  as 

fi{x,y,z,s,t)  =  x  -  hx(s,t)  =  0 
f2{x,y,z,s,t)  =  y-h2(s,t)  =  0 
f3(x,y,z,s,t)  =  z-h3(s,t)  =  0 

with  surface  point  x°  =  (x°,  y°,  z°,  a0,  t°).  The  gradients  of  the  /,  are 

V/i(x°,  y°,  z°,  s°,  t°)  =  (1,0,0,  hla,  hu) 

Vf2(x°,y0,z0,s°,t0)  =  (0, 1,0,  h2'3,h2't) 

V/3(x0,y0,z0,30,t°)  =  (0,0, 1,  h3iS,  h.3't) 

where  hU3  and  hut  are  partial  derivatives  of  h,  with  respect  to  s  and  t  respectively 
and  evaluated  at  (s°,  £°). 

The  linear  system 

C\  =  atihi'S  4-  a2h2'S  4-  a3h3'S  =  0 
C5  =  otihij  4-  ct2h2't  4*  ot3h3't  =  0 

has  a  solution 

(di,  q2,  d3)  =  (h2sh3(  —  h2th3a ,  huh3s  —  h\3h3t,  h\3h2t  —  hith23) 

which  is  {h\s,  h2s,  h33)  x  (hu,  h2t,  h3t).  Thus,  the  corresponding  C  as  defined  in  the 
algorithm  is  (di,  d2,  d3, 0, 0)  and  (dt,  d2,  q3)  is  the  normal  vector  of  the  surface  at 
(x°.y°,z0). 

The  D F(x°)  h  =  0  has  the  general  solution 

h23,  n3j,  —1.0)  4-  A?(  h\t ?  ^2*,  h3t,  0,  —  1 ) 

Thus  the  basis  vectors  of  the  affine  tangent  space  are 


hi  =  (hl3Ji23,h33, -1.0)  and  h2  =  (hu,  h2t,  h3t,  0. -1 ) 
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and  their  projections  into  (xj,  x2,  X3)-space 

hi  =  h.2S,  h,3a)  and  h2  =  (/iif, /i2*, /i3t) 

form  the  basis  of  the  corresponding  tangent  space  in  (x,y,  x)-space.  ■ 

Example  3.2  Consider  the  implicit  surface  g(x,y,z)  =  0  and  its  offset  given 
in  Example  1.1.  £>F(x°)  of  (1.4)  at  surface  point  x°  =  (x°,  y°,  :°,  u°,  v°,  w°)  = 
(0,0, 4,  0,0, 2)  is 

0  0  4  0  0  -4 

0  0  0  0  0  4 

2  0  0  -10  0  0 

0  2  0  0  -4  0 

The  linear  system  with  ai,  a2,  <*3,  a4  as  unknowns  is 

C4  =  -10q3  =  0 

Cs  =  — 4a4  =  0 

Cs  =  4a2  —  4qi  =  0 

and  has  a  solution  (dj,  a2,  a3,  a4)  =  (l,  1,0,0).  Hence  the  corresponding  C  in  the 
algorithm  is  (0,  0, 4, 0, 0,0)  in  which  (0, 0,  4)  is  the  normal  vector  at  (x°,  y°,  r°). 
The  general  solution  of  DF(x°)  h  =  0  is  Ajl^  +  A2h2  where 

h,  =  (1,0,  0.1/5, 0.0) 

and 

h2  =  (0,2, 0,0, 1,0) 

serve  as  the  basis  of  the  affine  tangent  space.  By  natural  projection, 

hx  -  (1,0,0) 

and 

h2  =  (0,2.0) 


are  the  basis  vectors  of  the  corresponding  tangent  space  in  (x.  y,  c  )-space. 
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3.1.3  Normal  Curvature 

For  a  2-surface  in  R3,  formulae  for  computing  the  normal  curvature  of  the  sur¬ 
face  along  a  tangent  direction  at  a  surface  point  can  be  found  in  most  differential 
geometry  books.  However,  we  wish  to  compute  the  normal  curvature  of  the  pro¬ 
jected  surface  5/,  in  (xi,  x2,  x3) -space,  from  a  given  2-surface  F(x)  =  0  in  Rn.  It 
is  clear  that  one  could  first  determine  the  local  explicitization  or  local  parameter¬ 
ization  (both  are  in  (xj,  x2,  x3)-space),  discussed  later  in  subsequent  sections,  and 
then  apply  the  standard  formulae.  In  this  section,  we  describe  a  method  that  di¬ 
rectly  computes  the  normal  curvature  of  the  surface  5/  from  the  high-dimensional 
surface  representation  F(x)  =  0,  without  constructing  a  local  approximant  first. 

3. 1.3.1  The  Normal  Curvature  of  Hypersurface  in  Rn 

Let  g(xi,X2, . . . ,  xn)  =0  be  a  hypersurface  in  Rn  and  let  Sg  represent  the  zero 
set  of  5  =  0.  For  p  €  Sg  and  v  €  Tp(Sa ),  and  a  normal  vector  field  N  on  Sg,  we 
defi.  e  the  linear  map  Lp  :  Tp(Sg)  —*  Tp(Sg)  as 

L„(v)  =  -VvN  =  -(:Vo  3)'(0) 

whc  e  3  :  I  — *  Sg  is  any  parametrized  curve  on  Sg  with  3(0)  =  p  and  3(0)  =  v.  and 
(.V  o  3)(<)  =  ,V(3(0)-  The  definition  makes  sense  since  the  directional  derivative 
VViV  is  tangent  to  Sg.  Lp  is  usually  called  the  Weingarten  map  or  shape  operator 
of  g  at  p,  see,  e.g.,  [50,  66 j.  The  generalization  of  Meusnier’s  theorem  to  high 
dimensions  states  that  3(0)  •  N{p)  is  invariant  for  all  parameterized  curves  3  on 
Sg  ith  3(0)  =  p  and  3(0)  =  v. 

Lemma  3.1  (Meusnier)  Let  3  be  a  parameterized  curve  on  the  hypersurface  g 
in  Rn  with  3(0)  =  p  and  3(0)  =  v.  Then 

Lp(  v)  •  v  =  3(0)  •  Mp) 
where  A’  is  a  normal  vector  field  on  Sg. 
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If  N  is  the  unit  normal  vector  field  on  S3  and  v  is  a  unit  tangent  vector  of  S3 
at  p,  the  normal  curvature  k3{v)  of  Sg  at  p  in  the  direction  v  is  defined  as 

Mv)  =  lp(v)  • v 

The  computation  of  Lp(v)  •  v  can  be  conducted  according  to  the  following 
lemma: 


Lemma  3.2  Let  N  =  Vg,  Hg  be  the  Hessian  matrix  of  g  at  p,  and  v  be  a  tangent 
vector  of  g  at  p.  Then 

LP(V)  •  v  =  -vTHgv 


Moreover,  when  v  is  a  unit  tangent 


Mv)  = 


II^(P)II 


'Lp(v)  •  v 


Proof: 


LP(v)  •  v  =  -VviV-v 


-VyVp  V 

_(Vv|2.,vv| L . vJt-Vv 

\  UX\  OX 2  OXjt ) 

4A*"- 

^  d*9  ,  x 

-  L,  3"“a — (p)vivi 
,,}=i  ox,dxj 

—vtH~v 


..£^<pk 


dx,dxn 


According  to  Meusnier's  Theorem,  when  v  is  a  unit  tangent  we  have 

i,(v)  =  ,J(0) '  ii^Hi 

where  3  is  a  parameterized  curve  on  S3  with  J(0)  =  p  and  J(0)  =  v.  Thus 


Mv)  = 


1 


il^7(P)ll 


Lp(v)  •  v 
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3. 1.3.2  The  Normal  Curvature  of  Projected  2-D  Surface  in  R3 

For  a  2-surface  Sf  in  Rn  and  its  projected  2-surface  5/  in  (xi,  i2>  Xaj-space,  we 
intend  to  compute  the  normal  curvature  of  Sj  in  a  given  tangent  direction  directly 
from  F(x)  =  0  without  computing  /. 

Let  x°  €  5/  and  v  6  I*o(S/)  be  the  projections  of  x°  €  Sf  and  v  €  Txo(Sf ), 
respectively,  into  (xj,  x2,  x3)-space.  In  addition,  let  N,  be  the  normal  vector  field 
on  fi  =  0.  Since  x°  €  Sp  and  v  is  a  tangent  to  Sf  at  x°,  x°  is  a  surface  point  of 
/,  =  0  and  v  is  a  tangent  vector  to  /,  =  0  at  x°  for  z  =  I,  2, . . . ,  n  —  2. 

Let  Qj, . . . ,  an_2  be  real  numbers  such  that  YStZi  a,:V,(x°)  =  (A,  B,  C,  0, - 0). 

where  Ax  +  By  +  Cz  +  D  =  0  is  the  tangent  plane  of  /  =  0  at  x°,  i.e.,  (A.  B.C)  is  the 
normal  vector  of  5/  at  x°.  In  the  following,  we  will  show  that  the  L^v)  •  v,  where 

L’xo(v)  =  — VvAr,,  fori  =  1 . n  —  2,  and  the  linear  combination  q,(Lx0(v)  • 

v)  play  an  important  role  in  computing  the  normal  curvature  of  f{xx,x 2,x3)  =  0. 

Consider  any  curve  3{t)  =  (di(0*  •  •  • «  3«(0)  on  Sp  with  /3(0)  =  x°  and  d(0)  = 
v.  Then  3  is  a  surface  curve  on  /,  =  0  as  well,  for  i  =  1. 2, ....  n  — 2.  By  Meusnier’s 
theorem  we  have 

I'xo(v)  •  v  =  <,3(0),  .V,(x°)>  (3.5) 

Hence 

£  a,(L^o(v)  •  v)  =  a,  <i(0).  V,(x°)> 

i=i  .=i 

=  £  <J(0).a,A,(x°)> 

1=1 

=  <J(0),£q.,V1(x°)> 

>=i 

=  <j(0).(A.S.C,0 . 0)> 

=  <(j,(o).  J2(0)J3(0)),(A.p.a> 

where  ( 3x(t),  J2(C-  h(t))  is  the  projected  curve  of  3  on  5/  with 

(■A(0).  J2(0).  Ji(0))  =  x°  and  l  .3,(0).  j2(0).  j3(0))  =  v 


(3.6) 


67 


and  (A,B,C)  is  a  normal  vector  of  Sf  at  x°.  Notice  that  (AXB,C)  and  v  may 
not  be  a  unit  normal  and  a  unit  tangent,  respectively,  to  5/  at  x°.  For  the  right 
hand  side  of  (3.6)  to  represent  the  normal  curvature  of  Sj,  both  the  normal  vector 
(A,  B ,  C)  of  S/  at  x°  and  (0i(O),  /?2(0),  /?3(0))  must  be  unit  vectors.  Thus,  when  0 
is  a  curve  such  that  (0i(O),  /?2(0),  /33(Q))  is  a  unit  vector,  according  to  Meusnier’s 
theorem,  the  right  hand  side  of  (3.6)  divided  by  the  length  of  ( A.B,C )  is  the 
normal  curvature  of  5/  at  x°  in  the  direction  (/?i(0),  ^(O), /?3(0)).  Hence,  the  left 
hand  side  of  (3.6)  divided  by  ||(.4,  B,C)\\  is  the  normal  curvature  of  Sj  at  x°  in 
the  direction  of  v,  provided  that  v  is  a  unit  tangent.  Notice  that  the  equality  in 
(3.5)  is  valid  for  a  non-unit  normal  vector  field  Nx  and  any  tangent  vector  v.  We 
summarize  this  fact  as  follows: 

Theorem  3.3  Let  N,  =  V/),  i  =  1, 2, ....  n  —  2  and  let  aj , . . . ,  an_2  be  numbers 
such  that  cttNt{x°)  =  (A,B,C,  0,...,0),  where  (A,  B,C)  is  the  normal  of  Sf 
at  x°.  Also  let  v  be  a  tangent  vector  of  Sf  at  x°  such  that  v  is  a  unit  tangent  to 
Sf  at  x°.  Then 

«xW£a'(LWv)v) 

is  the  normal  curvature  of  S j  at  x°  in  the  v  direction. 

The  theorem  and  its  proof  suggest  to  us  two  computation  schemes: 

Algorithm  3.2 

1.  Determine  v  such  that  v  is  unit. 

2.  Derive  3(t)  =  (0\(t), . . . ,  0n(t))  on  such  that  J(0)  =  x°  and  i(0)  =  v. 

3.  Compute  (3t( 0), ,ij2(0),  J3(0)). 

4.  Compute  (.4,  B.  C,  0 . 0)  =  EEV  a,.Y,(xu)  for  some  a x . an_2  and  nor¬ 

malize  (/ 4.  B.  C)  to  ( .4,  B.  Cl  with  unit  length. 

5.  Compute  <(it(0).  J2(0).  B.  < '  > 
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Algorithm  3.2  is  conceptually  very  simple,  however,  the  computation  of  the  3  on 
a  2-surface  Sf  is  nontrivial. 

Note  that,  as  shown  in  Lemma  3.2,  L^v)  •  v  =  — vr//j  v,  where  //,  is  the 
Hessian  of  the  hypersurface  /,  at  x°.  Moreover,  because  of  the  bilinearity  of  the 
form  vT Hi  v,  we  have 

£  a«(Lxo(v)  •  V)  =  Ct,(-VTH,  v)  =  -vT(]T  atHi)v 

.=i  i=i  i=i 

In  the  following,  we  state  an  algorithm  that  is  well  suited  to  computing  the  normal 
curvature  of  5/  at  x°  in  different  tangent  directions. 


Algorithm  3.3 


1.  Compute  Ni  =  Vfi,  for  2  =  1,2 . n  —  2. 

2.  Compute  a1,...,an-2  such  that  ET=i2  Oft;V,(x°)  =  (A,  £?,C.  0 . 0),  where 

(A,  B ,  C)  is  the  normal  of  5/  at  x°. 

3.  Compute  Ha  =  ociHi,  where  Hi  is  the  Hessian  of  /,  at  x°. 

4.  Adjust  v  such  that  v  is  a  unit  tangent  to  5/  at  x°. 

5.  The  normal  curvature  of  5/  at  x°  in  the  v  direction  is 


-1 


-vr//0  v 


11(^.5,011 

Let  L(q  be  the  shape  operator  of  the  closed  form  /(xl5 x2,x3)  =  0  at  x°.  Due 


to  Meusnier's  theorem,  from  (3.6)  we  obtain  a  formula  for  computing  L^0(v)  ■  v. 


Theorem  3.4  Let  N,  =  V/,,  i  =  1.  2, ....  n  —  2  and  let  a\, . . . .  an_ 2  be  numbers 

such  that  OiiV.fx0)  =  ( A,  B,  C,  0 . 0),  where  (A,  B.  C)  is  the  normal  of  S  j 

at  x°.  Also  let  v  be  a  tangent  vector  of  Sf  at  x°.  Then 

n—  2 

L^o( v)  •  v  =  Y,  a.iL^o(v)  ■  v) 

l=i 
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Example  3.3  Consider  the  implicit  surface  g(x,y,z)  =  0  and  its  offset  given 
in  Example  1.1.  Let  Ni  be  the  gradient  for  the  i-th  polynomial  /,,  i  =  1,2.  3, 4, 
x°  =  (i°,  y°,  z° ,  u°,  v°,w°)  =  (0,  0,4, 0,0, 2),  and  v  =  (0, 1,0,0, 1/2,0)  be  a  tangent 
vector  at  x°.  In  Example  3.2,  we  have  computed  (at,  a2,  a3,a4)  =  (1, 1,0,0)  such 


that 

(0, 0, 4, 0, 0, 0)  =  ^  aiNt(x°) 

1=1 

Now,  we  compute  Hq  =  a,//,  and  obtain 


2  0  0  -2  0  0 

0  2  0  0  -2  0 

0  0  2  0  0  -2 

-2  0  0  10  0  0 

0  -2  0  0  4  0 

0  0  -2  0  0  4 


Thus, 


-1 


-vtH0  v  =  -1/4 


11(0,0,4)11 

which  is  expected  to  be  the  normal  curvature  of  5/  at  (0,0,4)  in  the  direction 
of  (0,1,0).  Intersecting  the  projected  surface  with  the  plane  that  goes  through 
(0,0,4)  and  is  spanned  by  the  normal  (0,0,4)  and  tangent  (0.1,0),  we  find  that 
the  normal  section  is  a  circle  of  radius  4  and  has  normal  curvature  —1/4  at  (0. 0,  4) 
in  the  direction  of  (0, 1,0),  which  verifies  the  result.  ■ 


Example  3.4  Consider  the  Voronoi  surface  of  g  and  h  in  Example  1.2.  Let  .V,  be 
the  gradient  for  i-th  polynomial  /,,  i  =  1, 2. ....  8,  and  let 

x°  =  (i°,  y°,  z°,  u°.  v°.  w°,  u°,  u°,  w° .  r°)  =  (0, 0. 0. 0. 0, 1. 0. 0,  - 1, 1) 

and  v  =  (0, 1.0, 0. 1,0.0. 1/2, 0.0)  be  a  tangent  vector  at  x°.  We  find 


(ai,  a2 


as)  =  (1.  1.0.0. -1. -1.0.0) 
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3.2  Degree  Two  Implicit  Approximation 

Our  task  here  is  to  derive  a  degree  tv  o  implicit  approximant  of  5/  in  (i1?  x2,  x3)- 
space  for  the  given  F(x)  =  0  and  a  surface  point  x°.  We  require  that  the  approx¬ 
imant  has  the  same  normal  direction  as  Sr  at  x°  and  that  its  normal  curvatures 
agree  with  those  of  5/  in  all  tangent  directions  at  xu.  The  constraint  on  curvature 
agreement  in  all  directions  is  difficult  to  r<\i  use  straight  forwardly.  However,  using 
the  following  Three  Tangents  Theorem  1 52.  arc  .mie  to  tormulate  the  constraint 
easily. 
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Theorem  3.5  (Three  Tangents  Theorem)  Two  surfaces  with  common  normal 
direction  at  a  surface  point  have  the  same  normal  curvatures  in  every  tangent  di~ 
rection  if  and  only  if  their  normal  curvatures  agree  in  three  tangent  directions  of 
which  any  two  are  linearly  independent. 

Applying  this  theorem,  the  constraints  above  can  be  translated  into  a  system 
of  seven  equations,  some  linear  and  some  nonlinear.  The  nonlinear  equations  stem 
from  the  normalization  of  the  normal  vector  of  the  approximant,  in  the  normal  cur¬ 
vature  computation.  Using  Lemma  3.2  to  rephrase  the  Three  Tangents  Theorem, 
however,  results  in  a  system  of  linear  equations  which  enforce  the  constraints. 

Corollary  3.1  Two  surfaces  with  identical  normal  vector  at  a  surface  point  p  have 
the  same  normal  curvatures  in  every  tangent  direction  if  and  only  if  values  of  their 
Lp(v)-v  agree  in  three  tangent  directions  of  which  any  two  are  linearly  independent. 

Let  g2(xi,x2,x 3)  =  Hl+j+k=oaijkxix2x3  be  a  generic  degree  two  polynomial 
with  the  symbolic  coefficients  aq*.  Suppose  vx,  v2  and  v3  are  tangents  to  5/r  at  x° 
with  the  property  that  any  pair  of  their  projections,  vx,  v2  and  v3  respectively,  in 
(x\,  x2,  x3)-space  is  linearly  independent.  We  translate  the  constraints  to  a  linear 
system  as  follows: 

1.  r(x°)  =  o. 

2.  V<72(x°)  =  (A,  5,C),  where  (A,  B.C)  is  the  normal  vector  of  5/  at  x°  com¬ 
puted  using  Algorithm  3.1. 

3.  For  i  =  1, 2, 3, 

L^o(v.)  •  v,  =  L(0( v,)  •  v, 

where  L^0(v,)  •  v,  =  and  L(0(v,)  •  v,  is  computed  according  to 

Theorem  3.4. 

This  formulation  results  in  7  linear  equations  in  10  variables.  Since  g:(x\.  x2,  x3)  = 
0  and  yg~\ Ji.  x2,  x3)  =  0.  where  ^  0.  represent  the  same  surface.  x2,  x3) 
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has  9  coefficients  on  which  the  surface  depends.  Thus,  we  can  adjoin  to  the  linear 
system  the  additional  equation 

atJk  -1=0 

for  some  0  <  i  +  j  +  k  <  2,  and  obtain  a  8  x  10  linear  system  from  which  the 
coefficients  of  g 2  can  be  computed  with  two  degrees  of  freedom.  We  give  the 
algorithm  as  follows: 

Algorithm  3.4 

1.  Form  a  generic  polynomial  ^(x, ,  ,r2,  x3)  =  T.‘i+j+k=oaijkx\x}2x$,  where  the 
a,,*,  are  symbolic  coefficients. 

2.  Compute  N,  =  V/,,  for  i  =  1, 2. ....  n  —  2. 

3.  Compute  Q1,...,Qn_2  such  that  £?=i2  atN,(x°)  =  (A,B,C.O . 0),  where 

(.4.  B,  C)  is  the  normal  of  5/  at  x°. 

4.  Compute  Hq  —  £?_T2  a,#,,  where  H,  is  the  Hessian  of  /,  at  x°. 

5.  Derive  three  tangent  vectors  v,,  v2  and  v3  to  Sf  at  x°  such  that  any  pair  of 
their  projections.  Vi,v2  and  v3  respectively,  in  (xt, x2, x3)-space  is  linearly 
independent. 

6.  Let  V<72(x°)  =  gl^x0),  g^{x0)). 

7.  Form  the  following  linear  system: 

a,jk  —  1  =0.  for  some  0  <  i  4-  j  +  k  <  2 

<?2(x°)  -  0 

/72,(x°)  -  A  =  0 

gl7(x°)  -  B  =  0 

1721(x0)-C  =  0 

vf H3vt  -  vf  H0  v,  = 


0 
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vj Hgv2  -  ^7  Ho  V2  =  0 
v^Hgv3  -  V3  H0  v3  =  0 

8.  Solve  the  above  system  for  the  coefficients  of  g 2. 

Example  3.5  Consider  the  implicit  surface  g(x,y,z)  =  0  and  its  offset  given 
in  Example  1.1.  Let  x°  =  (x°,  y°,  z°,  u°,  v°,  w°)  =  (0,  0,  4, 0,  0, 2),  and  let  Vj  = 
(1,2,0, 1/5, 1,0),  v2  =  (1,1,  0,1/5, 1/2,0),  v3  =  (2. 3,  0, 2/5, 3/2, 0)  be  three  tan¬ 
gents  at  x°.  Note  that  any  pair  of  their  projections  is  linearly  independent.  Let 
g2(xi,x2,x 3)  =  H,2+j+*r=o a,jkx\xJ2x^  be  a  generic  degree  two  polynomial  with  the 
symbolic  coefficients  a We  form  the  following  linear  system: 

16aoo2  4-  4  a  ooi  4-  a  000  —  0 
4aioi  +  aioo  =  0 

^aon  +  a0io  =  0 

8  a  002  +  a  001  —  4  =  0 
— 2a2oo  ~  -Quo  —  80020  —  -ano  4*  28/5  =  0 

— 2a20o  —  2au0  —  2ao2o  4-  13/5  =  0 

—8a 2oo  —  6ano  —  6Qno  ~  lSao2o  +  n  /5  =  0 

When  aooo  —  1  =  0  is  added  to  the  system,  for  the  unknown  coefficients 


(aooOi  QlOOi  flOlOi  QOOI 1  a200*  ^HOi  OlOli  ^020i  ami ,  0002 ) 
we  obtain  the  general  solution 


( 1,  -4a,  -4 3.  -9/2,  4/5,  0.  a,  1/2,  3 , 17/16) 


With  a  =  3  =  1,  we  have 


g2(x,y.z)  =  1  -  4x  -  4 7 


9  4  2  1  „ 

4-  4r:4  -y~  +  yz  + 

1  •)  2 
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3.3  Local  Explicit  Approximation 

The  implicit  function  theorem  ensures  that  system  (3.1)  determines  n  —  2  com¬ 
ponents  of  x  as  functions  of  the  remaining  2  components  in  a  neighborhood  of  any 
surface  point  x°  at  which  DF(x°)  has  rank  n  —  2.  Since  D F(x°)  has  maximal 
rank  n  —  2,  some  set  of  n  —  2  columns  of  its  matrix  is  linearly  independent.  In 
the  following,  without  loss  of  generality,  we  assume  that  the  last  n  -  2  columns 
are  linearly  independent,  and  hence  the  determinant  of  the  Jacobian  matrix  of 
/i,  /2, . . . ,  fn-2  with  respect  to  x3,  x4, . . . ,  xn,  denoted  as  j  JF(x°)|,  is  nonzero.  The 
implicit  function  theorem  guarantees  that  there  exists  a  neighborhood  V  of  x°.  an 
open  set  U  C  E2  containing  (xj,x£),  and  a  mapping  $  =  (03,  rpA, ....  wn)  denned 
on  U  such  that  |JF(,x)|  ^  0  for  all  x  €  V,  and  for  1  <  i  <  n  —  2. 


/,(X1,X2,V3(X1,X2),....C>„(X1,X2))  =0  (3.7) 

for  all  ( x  i ,  x2 )  €  U.  Let  tp{  be  defined  as 

0,(xi,x2)=  b\'kx\A  (3.8) 

j+fc>  o 

and  /.(i!,^)  denote  /.(x,,  x2,  v3(xx,  x2), - iyn(xt,  x2)). 

It  is  clear  that,  in  (xt,  x2,  x3)-space,  x3  —  v3(xi .  x2)  represents  the  local  explicit 
approximation  of  the  projected  surface  (3.2)  at  (x°,x^,  x3),  and  might  be  called  a 
local  exphcitization  of  the  surface.  The  local  explicitization  of  an  implicit  surface 
g{x i,x2,x3)  =  0  has  been  considered  in  [48]. 

The  unknown  coefficients  of  can  be  calculated  from  the  chain  rule  with  a 

linear  system  solver.  By  requiring  that  the  k- th  derivatives  of  /,,  i  =  1 . 2, _ n  —  2. 

witn  respect  to  xt  and  x2  are  identically  zero,  a  linear  system  with  unknown 
coefficients  of  degree  k  terms  is  obtained.  It  is  structurally  very  simple.  However, 
the  direct  application  of  the  chain  rule  results  in  a  formula  which  is  algebraically 
tedious.  In  the  following,  we  develop  a  recursive  formulation  which  presents  the 
computation  in  a  more  suitable  manner. 


75 


When  we  assume  that  x°  =  (0,0 . 0)  is  the  origin,  the  constant  terms  of 

/,,  doo  o,  and  the  constant  term  of  t/>„  6™,  are  both  zero  and  the  partials  -^4— r 

dr^dx^ 

evaluated  at  (0, 0)  and  divided  by  z!  and  ;'!  are  the  coefficients  of  xjx%  in  /,(xi,  x2). 
With  this  property,  a  recursive  formula  is  derived  as  follows. 

Let  matrix  A  be  the  Jacobian  matrix  of  ft, , . . . ,  /„_2  with  respect  to  x3, , . . . ,  x„ 
evaluated  at  x°,  i.e. 


A  =  J  F(x°)  = 


a(1) 

a0010...0  a00...01 


(n-2)  (n-2) 

“0010... 0  “00...01 


We  also  define  (d,(j))kt  as  in  [4S],  which  is  the  coefficient  ofx^Xj  in  (t/»i(x1,x2))J, 
that  is, 

(V,(x1,x2))J  =  £  (di(j))ux\xl7 

k+l>j 

Analogous  to  the  recursive  derivation  for  {a(k))i}  in  Section  2.2.1,  the  {d,(j))ki  can 
be  recursively  defined  as 


1  if  j  =  k  =  /  =  0 

......  0  if  j  =  0  and  k  +  l  >  1 

(di(j))ki  =  . 

if  j  =  1  and  k  +  l  >  1 
,  —  l<p+q<k+l-l  ( di(j  1  ))pq^(k-p)(i~q)  if  J1  ^2  and  k  +  l  >  j 

The  formula  for  computing  6^,  where  j  +  k  >  1  and  z  =  3.  4, _ n,  is: 

For  j  +  k  =  1 


’’’Ji  :  ,I')io...o 
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For  j  +  k  >  2 


where 


.<•')  _  „(0 


=  ajkO...O 

j+k—l  kt+l  3 

+  E  E  E 

+J2=0  *3+  ■  +  *««  j-ii  J3=0 

+  s  '  — 


...j„  (<^3( J3)  ) ^3/3  ‘  '  '  (dn(jn))knt 


J3+-+Jn^2-Ji  -J2 


Example  3.6  Consider  the  surface  g  =  0  and  its  offset  given  in  Example  1.1  and 
the  surface  point  x°  =  (x°,  y°,  20,  u°,  v° ,  u;0)  =  (0,  0, 4, 0,  0, 2).  First  of  all.  we  have 
to  translate  x°  to  the  origin,  and  as  the  result,  system  (1.4)  becomes 


z2  —  2  wz  +  y2  —  2vy  +  x2  —  2ux  +  w2  +  v2  +  u2  +  4a  —  4ih  =  0 

w2  +  v2  +  4u2  -f  4u;  =  0 
— 4ui  4-  wx  +  3uu>  4-  2i  —  lOu  =  0 
—  vz+wy  r  2y  —  4u  =  0 

which  has  the  Jacobian  matrix  A  with  respect  to  z,u,v,w  and  evaluated  at 

(0,0,0.  0.0,0), 

4  0  0  -4 

0  0  0  4 

A  = 

0  -10  0  0 

0  0-40 

Solving  two  linear  systems,  we  have  coefficients  for  linear  terms 


f-(3)  ,(•*)  ,(5)  ,(6K 
'  °10  1  °10  >  °10  '  °10  / 

= 

(0, 1/5, 0,0) 

,.(3)  ,(■*)  ,<5)  ,(SK 
(  ®01  !  1  ’  °01  ) 

= 

0 

cT 

r  *4 

O 

O 

we  have 

1  J>>  J2)  „n)  n, 

(  ~r20  1  —  r20  1  ~ r20  •  r  JO  ) 

=  (-1. 0,0.0) 

_,(2)  _rn>  _r 

'  rl  1  ’  rU  '  rn  •  r 

1 1  ’  * 

=  (0,0.0. 0) 

1 

*o  — * 

1 

"1 

1 

rsJ  >.«. 

1 

h 

■  ,2  1 

=  !  -1.0..  0.0) 

I  i 


and  have  coefficients 

(4o,4o  ,4o\©  =  (  —  1/4, 0,0,0) 

(^n^n^u^u)  =  (0,0, 0,0) 

(C^otC,^)  =  (  —  1/4, 0,0,0) 

Thus,  z  =  —  (x2  +  y2)/4  and  after  the  reverse  translation  2  =  4  —  (x2  +  y2)/4  is  the 
degree  2  local  explicitization  of  the  surface  at  (0, 0,  4, 0, 0,  2).  ■ 

3.4  Local  Parametric  Approximation 

In  [14],  a  Taylor  approximant  of  the  intersection  curve  of  two  implicit  surfaces 
is  constructed  which  serves  as  the  local  parametrization  of  the  intersection  curve. 
We  generalize  this  method  to  our  problem  domain. 

For  a  given  system  of  equations  (3.1)  and  a  point  x°  on  it,  we  seek  a  paramet¬ 
rically  described  solution 

<&(u,v)  =  (<Mu,v),<Mu,v),...,0n(u, v))  with  4>(0,0)  =  x°  (3.9) 

in  the  neighborhood  V  of  x°.  Solution  (3.9)  is  basically  a  local  parametrization  of 
the  surface  represented  by  (3.1)  at  point  x°.  The  three  coordinate  functions 

Xi  =  o1(u,v),  x2  =  <z>?(u,  v),  x3  =  o3(u,v) 

define  a  local  parametrization  of  the  projected  surface  (3.2)  in  (xt,  x2,  x3)-space. 

When  the  hvpersurfaces  /,  intersect  transversally  and  are  not  singular  in  the 
vicinity  of  x°  such  a  parametrically  defined  solution  exists.  In  the  neighborhood 
of  x°.  V,  the  surfaces  are  hence  defined  as  the  solution  of 

f,(u.  v)  =  /,(<3>i(u,  t’).  02(u,v) . on(u.  V))  =  0.  i  —  1,2, - n  -  2. 

The  Taylor  expansion  of  /,(u,  v)  in  power  of  u  and  v  is 


/.( u-  L') 


/,(0.0)  4-  u 


(1L 

du 


r)f. 

'  <)v 


<r_cH 

2  ()u: 


OvcJv 


±d±L 

2  dvl 


(3.10) 
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=  /,(*(0,0))  +  u 


I v 

v£;  dxi  du ) 


+  V 


( V'  c^/i  d<t>]  \ 

\h  dx>  dv  ) 


d2f,  dd>}\  d^]  df^PokX  , 
dx:dxk  du  J  du  dik  du 2  J 

=  /,(*(0.0))  +  UV/,(x°)  •  *.(0,0)  +  «v/,(x°)  •  *„(0,0) 
+  y(v/,(x0)  •  *uu(0,0)  4-  *.(0,0)  •  ///.  •  *.(0,0)] 


+  uv[V/,(x°)  •  *..(0, 0)  +  *„(0, 0)  •  Hft  •  *u(0, 0)] 
+  y[V/,(x°)  •  *uu(0, 0)  +  *u{0, 0)  •  Hfi  ■  *„(0,  0)] 


+  ••  • 


(3.11) 


where  Hft,  the  Hessian  of  /,,  and  all  the  derivatives  of  /,  are  evaluated  at  x°.  and 
all  the  derivatives  of  /,  and  <t>,  are  evaluated  at  (0,0). 

In  equation  (3.10).  the  coefficient  of  each  monomial  uJvk  must  vanish,  so  by 
requiring  the  coefficient  of  uJvk  be  zero  for  1  <  j  +  k  <  m.  a  truncated  local 
parametrization  of  degree  m  can  be  derived.  This  amounts  to  solving  the  following 
series  of  linear  systems 

V/,(x°)  •  *uJv*(0,0)  =  B^,  i  =  1.2 . n-  2  (3.12) 

where  1  <  j  +  k  <  m.  *uJu*(0.0)  =  ( jj~r(0. 0)) ^  ^  .  and  the  are 
expressions  of  partial  derivatives  of  /,  and  lower  degree  partial  derivatives  of  o,, 
for  example. 


d(  *) 
"10 

=  0 

(3.13) 

o(>) 

UQ\ 

=  0 

(3.14) 

D<  *) 
&20 

=  -*,,(0.0)  ■ 

//,. 

*.,(0.0) 

(3.15) 

B\‘! 

=  -*„(0.0> 

//,. 

*.,(0.0) 

(3.16) 

ni  <) 
&02 

=  -*,,(0.0 

//.. 

*  tO.Oi 

(3.17) 

Note  that  the  B\'±  can  be  recursively  define.; 

>  *  he  origin  1 0. .  . 

.  .  0).  in  analogy 

to  Section  3.3. 


The  system  (3.12)  is  a  (n  —  2)  x  n  system.  Since  we  assume  that  Z?F(x°)  has 
rank  n  —  2,  the  solution  of  this  system  can  be  written  as  the  linear  combination 
of  ti,  t2,  V/i(x°), . . . ,  V/n_2(x°),  where  tt  and  t2  form  a  unit  vector  base  of  the 
tangent  space  at  x°,  that  is, 

0)  =  +  ag’t,  +  +  ■■■  +  (3.18) 

Substituting  into  (3.12)  yields 

x°)  •  V/,(x°)  =  Bj?,  .  =  1,2 . n-2.  (3.19) 

/=  1 

from  which  the  unique  solution,  . . . ,  can  be  derived.  Hence,  the  ex¬ 

pression  (3.18)  with  0$  and  ct^  arbitrarily  selected  is  the  general  solution  of  the 
system  (3.12).  The  suitable  values  for  and  q^'  are  chosen  as  follows: 

1.  For  systems  (3.13)  and  (3.14),  we  have 


We  assign  =  1  and  a**  =  0  for  (3.13)  and  aj*  =  0  and  =  1 
for  (3.14)  so  that  the  isoparametric  curves  $((u,0))  and  $((0,  v))  intersect 
transversally  at  point  x°. 

2.  For  systems  (3.15),  (3.16)  and  (3.17),  we  assign  =  a*,^  =  0. 

Note  that  the  three  coordinate  functions 


Jj  =  <Z>i(u,  v),  X 2  =  02(u,l’),  X3  =  03(u,l>) 

obtained  in  this  way  define  a  local  parametrization  of  the  surface  (3.2)  in  ( x,,  x2.  x3)- 
space. 

To  compute  the  solution  of  system  (3.12)  with  a  numerically  stable  method, 
we  consider  singular  value  decomposition,  see  e.g.,  [29].  The  matrix  (DF(  x°))r  is 


factored  as 


so 


where  U  =  [uj, . . . ,  u„]  €  Rnxn  and  V  =  [vj, . . . ,  vn_2]  €  R(n  2)x(n_2)  are  orthog¬ 
onal  matrices,  and  E  G  Rnx(n_2)  is  a  diagonal  matrix.  With  direct  substitution  of 
the  factorization,  system  (3.12)  becomes 

VETUT$uJvk{0, 0)  =  B 

where  B  =  [B^ , . . . ,  B^~2)]T.  Its  solution  can  be  generally  written  as 

$uJv*(0,0)  =  7iU!  +  72U2  H - H  7nU„ 

and  because  the  rank  of  the  differential  matrix  is  n— 2,  E^  ^  0  for  j  =  1,2, ....  n  —  2, 
and  thus  the  7lt . . . ,  7„-2  are  uniquely  determined  by 

7,  =  (vf  B)l  E„ 

and  7n_i  and  7„  are  arbitrary.  Here  Ui,  u2, - un_2  span  the  range  of  (DF(x°))r. 

hence  span  the  same  space  as  the  gradients.  Note  also  that  the  null  space  of  Z?F(x°) 
spanned  by  un_i  and  u„.  Thus  we  obtain  $u(0,0)  =  un_i  and  $„(0,0)  =  un.  For 
$u2(0,0),$uu(0,0),  and  $^(0,0),  we  let  7„_l  =  7n  =  0. 

Example  3.7  Consider  the  surface  g  =  0  and  its  offset  given  in  Example  1.1 
and  the  surface  point  x°  =  (x°,  y° ,  c°,  u°,  v°,  w° )  =  (0  0,  4, 0,  0, 2).  It  is  clear  that 
$(0.0)  =  1 0.0.  4, 0,0, 2),  $,(0,0)  =  tj  and  $,(0,0)  =  t2,  where  tt  =  (5/\/26. 0.0. 1/v '2(5.0.  O' 
and  t2  =  (0, 2/\/5,  0, 0,  l/\/5, 0).  By  applying  formula  (3.19),  we  have 

(4o,4o\4o\4o)  =  (-5/52,-3/26,0,0) 

(4n\4?  ,4?  >4?)  =  (0.0.0.0) 

(4J\  $£,  4?,  4?)  =  (-1/20.-3/40.0.0) 

Thus,  from  formula  (3. IS),  it  follows  that 

$,j(0,0)  =  (0,0. -5/13.0,0. -1/13) 

$,*(0.0)  =  (0.0.0. 0.0.0) 

$fI(0.0)  =  (0.0. -1/5.0.0.-1/10) 


SI 


Consequently,  we  have 


<MM)  = 

<MM)  =  75  * 

<MM)  =  4-^2~i^2 


as  the  local  parametric  approximant  of  degree  2  in  (xx,  i2,  x3)-space. 
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4.  PIECEWISE  APPROXIMATIONS  OF  2-D  SURFACES 

Implicit  surfaces  have  recently  become  more  important  in  CAGD  and  solid 
modeling.  In  part,  implicit  surfaces  have  specific  advantages  over  the  traditional 
parametric  surfaces.  For  example,  many  complex  objects  can  be  modeled  more 
easily  using  implicit  surfaces,  and  certain  geometric  operations,  e.g.,  membership 
classification  problem,  can  be  handled  straightforwardly  when  implicit  surface  rep¬ 
resentations  are  available.  Moreover,  using  implicit  surfaces,  a  number  of  sophis¬ 
ticated  techniques  have  been  proposed,  e.g.,  the  substitution  blending  surfaces  of 
[34,  35,' 36,  37]. 

As  the  role  of  implicit  surfaces  increases  in  importance  in  solid  modeling  and 
CAGD,  rendering  implicit  surfaces  efficiently  becomes  a  crucial  support  in  surface 
design.  In  computer  graphics,  many  surface  rendering  algorithms  rely  on  piece- 
wise  linear  approximations  (PLAs)  of  a  surface  since  a  PLA  allows  one  to  take 
advantage  of  hardware  capabilities  and  reduces  the  cost  of  expensive  ray  casting 
in  the  rendering  process.  However,  while  the  PLA  of  parametric  surfaces  has  been 
extensively  studied  and  utilized  as  a  tool  for  the  evaluation  of  surface  intersections 
[10,  16,  38,  41,  42]  and  for  rendering,  it  seems  that  much  less  attention  has  been 
paid  in  the  literature  to  the  PLA  of  implicitly  defined  surfaces.  Recently,  Bloomen- 
thal  [17]  has  proposed  an  algorithm  for  computing  the  PLA  of  an  implicit  surface 
based  on  space  subdivision  using  octrees.  In  [9,  7,  8,  6],  a  simplicial  continuation 
algorithm  is  presented  for  obtaining  a  PLA  to  a  component  of  an  implicit  surface. 
Both  methods  fundamentally  rely  on  vertex  evaluation.  Rheinboldt  [56.  57]  has 
presented  an  algorithm  that  maps  a  triangulation  of  Rp  to  a  p- manifold,  where 
p  >  1,  and  hence  induces  a  triangulation  on  the  p- manifold. 
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Offset  and  Voronoi  surfaces  of  implicit  surfaces  can  be  formulated  as  the  projec¬ 
tion  to  R3  of  2-surfaces  that  are  implicitly  defined  in  six  and  ten  dimensional  spaces 
respectively.  Algorithms  using  space  subdivision  or  simplicial  continuation  can  be 
generalized  to  compute  the  PLA  of  the  2-surface  in  high-dimensional  space.  Thus, 
to  compute  the  PLA  of  the  offset  of  an  implicit  surface,  we  could  compute  the  PLA 
of  the  2-surface  in  6-space,  and  then  project  it  into  3-D  subspaces.  However,  as  the 
formulation  dimension  increases,  the  complexity  of  computing  the  PLA  increases 
exponentially.  To  reduce  the  complexity,  we  propose  an  algorithm  that  determines 
the  PLA  of  the  projection  to  R?  of  the  2-surface  defined  in  high-dimensional  space, 
with  all  major  computations  performed  in  3-space. 

In  Section  1  we  briefly  describe  a  method  based  on  the  simplicial  continuation 
method  due  to  AUgower  and  Gnutzmann  [7],  and  explain  some  difficulties  when  ap¬ 
plied  to  offset,  Voronoi  and  blending  surfaces.  In  Section  2  we  sketch  the  proposed 
algorithm  in  pseudo  code,  and  then  in  Section  3  we  describe  the  computations  in 
detail. 


4.1  The  PLA  of  Offset,  Voronoi  and  Blending  Surfaces 

As  described  in  Section  1.1.3,  the  offset  of  an  implicit  surface  can  be  viewed  as 
the  projection  to  R3  of  a  2-surface  in  R6,  and  the  Voronoi  surface  of  two  implicit 
surfaces  as  the  projection  to  R3  of  a  2-surface  in  R10  as  shown  in  (1.3)  and  (1.5) 
respectively.  The  2-surfaces  are  represented  by  the  following  system  of  n  —  2 
equations  in  n  variables, 

/l(il,z2,....xn)  =  0 
h(X\,Xl . -Tn )  =  0 

(4.1) 

In- 2(xli  •  •  •  ■ -^n)  =  0 

abbreviated  in  matrix  form  as  F(x)  =  0.  We  suppose  that  the  target  surface  we  are 
interested  in  is  in  (xi,  x2,  x3) -space.  A  closed- form  /(xj.  x2,  x 3)  =  0  can  be  obtained 
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in  principle  by  eliminating  the  last  n  —  3  variables,  but  this  is  often  not  possible 
in  practice.  It  is  also  impractical  to  compute  the  PLA  of  the  surface  in  n-space. 
using  generalizations  of  space  decomposition  or  simplicial  continuation,  because 
the  complexity  of  the  computation  increases  exponentially  with  the  dimension  of 
the  space.  One  way  to  reduce  the  complexity  might  be  by  computing  the  PLAs 
of  the  desired  surface,  e.g.,  offset,  and  its  basis  surface  in  parallel  in  such  a  ./ay 
that  the  major  computation  is  performed  in  the  orthogonal  3-D  subspaces.  We 
have  looked  into  such  an  approach  for  computing  the  PLA  of  an  offset  surface.  Let 
system  (1.3)  define  the  offset  in  R1 2 3 4 * 6.  For  a  pair  of  corresponding  points  and  x2 
on  the  offset  and  its  basis  surface  respectively,  we  determine  two  3-simplices.  one 
containing  xt  in  (x,  y,  x)-space  and  the  other  containing  x2  in  (u,  v,  ty)-space.  Then 
two  sequences  of  transversal  simplices  are  constructed  in  parallel.  T wo  sequences 
of  simplices  are  constructed  that  are  coordinated  by  the  following  computations 
and  considerations: 

Let  (7 1  and  a2  be  two  simplices  that  we  are  currently  processing,  Xj  6  ax  be  a 
point  on  the  offset  surface  and  x2  €  r2  be  the  footpoint  of  Xi  on  the  basis  surface. 

1.  Determine  the  simplex  a  C  Rb  such  that  (Xi,x2)  6  o  and  a  projects  to  <tx 
and  <r2. 

2.  Compute  the  affine  map  H  that  interpolates  F(x)  at  all  the  vertices  of  a. 

3.  Project  the  affine  map  H  to  Ha  and  H2  in  the  two  orthogonal  subspaces. 

4.  Compute  the  intersections  of  H,  with  edges  of  cr,,  for  t  =  1,2. 

•5.  Pair  the  intersections  with  ax  and  <72  to  obtain  points  in  6-space  and  refine 
them  iteratively  back  to  zeros  of  F(x)  =  0.  ^he  points  will  be  the  new 
(Xi,X2). 

6.  Based  on  Xi  and  x2,  we  determine  how  ro  pivot,  the  current  pair  of  simplices. 
The  resulting  simplices  are  the  : : » •  \ v  -r,  and  <r2. 
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7.  Go  back  to  Step  1. 

We  have  investigated  such  an  approach  and  found  it  unattractive  for  the  following 
reasons: 

•  For  given  a\  containing  Xi  in  (x,  y,  x)-space  and  c2  containing  x2  in  (u,  u,  w)- 
space,  there  are  4!4!  possible  simplices  a  in  6-space  that  project  to  the  same 
pair  <r i  and  cr2.  Hence  Step  1  requires  considerable  computation  for  generat¬ 
ing  a  and  checking  if  (xt,  x2)  €  o. 

•  For  two  transversal  simplices  <ri  and  <t2,  the  simplex  a  fo  ,nd  in  Step  1  might 
not  be  transversal.  Suppose  (bi,b2)  is  the  pair  of  simplices  in  the  previous 
step,  corresponding  to  the  simplex  b  in  6-space.  That  is.  (j\  and  cr2  are 
the  result  of  Divoting  b\  and  b2  in  3-space,  respectively,  with  respect  o 
shared  transversal  faces.  Note  that  the  pivoting  of  bx  and  b2  corresponds 
to  a  sequence  of  pivotings  of  b  which  is  not  necessarily  done  with  respect  to 
transversal  faces.  Therefore,  a  might  be  not  transversal.  Now  (xi,x2)  €  cr. 
and  so  it  is  easily  seen  that  in  such  a  case  the  surface  penetrates  a  although 
all  vertices  of  a  have  the  same  signs  on  vertex  evaluation.  In  this  case,  the 
computed  affine  map  H  is  useless  for  continuing  the  computation. 

•  The  two  projected  affine  maps  Hi  and  H2  are  generally  not  the  affine  maps 
of  (Ti  and  cr2  respectively.  Hence  the  projected  affine  maps  for  adjacent  3- 
simplices  are  often  not  continuous  at  the  shared  face. 

•  The  intersections  of  H,  with  a,,  for  i  —  1.2,  are  often  not  in  correspondence, 
and  then  it  is  not  possible  to  define  a  meaningful  point  pairing. 

•  When  applied  to  offsets  of  parametric  surfaces.  Steps  4.  5.  and  6  require 
substantial  modifications. 

Using  3-simplices  is  attractive  because  with  it  vertex  evaluation  is  unambiguous 
and  there  is  a  simple  pivoting  scheme.  However,  adaptive  subdivision  resulting 
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again  in  3-simplices  is  more  complicated.  So,  if  vertex  evaluation  is  not  necessary, 
using  cubes  to  subdivide  space  might  be  a  better  choice  since  adaptive  subdivision 
is  very  simple. 

The  parametric  approximant  4>(s,  t),  discussed  in  Section  3.4  could  be  used 
in  place  of  an  affine  map.  However,  the  coordination  of  the  sequences  of  cubes 
is  complicated,  especially  when  computing  Voronoi  and  blending  surfaces  where 
more  than  two  sequences  of  cubes  have  to  be  put  into  correspondence. 

4.2  Proposed  Approach 

Instead  of  computing  a  PLA  of  the  2-surface  in  n-space,  or  its  projection  into 
several  subspaces,  we  will  compute  the  PLA  of  one  projection  only.  In  the  case 
of  offset  surfaces,  this  means  that  we  approximate  the  offset  in  3-space  without 
explicitly  approximating  the  basis  surface  as  well,  or  the  2-surface  in  6-space  corre¬ 
sponding  to  both  the  offset  and  its  basis  surface.  We  will  use  a  method  similar  to 
simplicial  continuation  in  [8],  but  based  on  local  parametric  approximation.  The 
approach  is  well  suited  to  those  2-surfaces  that  have  been  defined  in  Rn,  but  whose 
projection  to  3-space  is  ultimately  of  interest. 

Let  F(x)  =  0  be  the  surface  definition  in  Rn,  5/  its  projection  into  (xj,  x2,  x3)- 
space.  Given  a  regular  surface  point  x  =  (xj,x2)  €  Rn,  where  Xi  €  R3  and  x2  € 
Rn~3,  and  a  cube  in  R3  containing  Xi,  a  sequence  of  consecutive  cubes  intersecting 
the  surface  5/  is  generated  in  a  fashion  that  spirals  outward  from  Xj.  For  each 
cube,  a  degree  two  parametric  approximant  $(s,  t)  =  (d>i(s,  f),  <£2(s,  t), . . . .  <?„(s,  <)) 
of  F(x)  =  0  is  derived,  and  $i(s,t)  =  t),  <f>2(s,  t),  <£3(s,  t))  serves  as  the 

surface  approximation  of  the  target  surface  in  (xt,  x2,  x3)-space.  The  parametric 
approximant  <$(.s,  t)  plays  an  essential  role  in  this  approach.  First  of  all.  the 
intersections  of  $i(s,  t)  with  the  edges  of  the  cube  are  used  as  estimates  that  are 
refined  to  surface  intersections  along  edges  of  the  cube.  Secondly,  the  parametric 
approximant  provides  a  way  of  recovering  a  point  in  Rn  while  doing  computations 
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in  R3.  Thirdly,  when  extending  the  surface  approximation  into  a  neighboring 
cube,  $(3,  t)  provides  a  mechanism  for  obtaining  an  estimate  to  a  surface  point 
x  =  (xj,  x2)  €  Rn  such  that  Xi  lies  in  the  neighboring  cube.  With  this  parametric 
approximant,  we  are  able  to  march  on  the  surface  in  (xt,  x2,  x3)-space  although 
the  surface  is  formally  defined  in  Rn,  n  >  3.  This  approach  is  of  course  also 
applicable  to  offsets,  Voronoi  surfaces  and  blends  where  some  of  the  basis  surfaces 
are  in  parametric  form.  For  the  computation  of  parametric  approximation,  see 
Section  3.4.  We  give  a  high-level  description  of  the  approximation  algorithm. 

Algorithm  4.1 

Input 
F(x)  =  0 

x  =  (xi,x2)  €  Rn:  a  regular  surface  point  where  Xj  6  (xt, x2, x3)-space 
B:  a  cube  with  edge  size  <5  containing  X! 

D:  a  compact  domain  (aj  :  a2,6j  :  62,cl  :  c2)  in  (x1?  x2,  x3)-space 
Output 

L:  List  of  transversal  cubes  and  their  surface  intersections  where 
the  cubes  are  contained  in  D  or  intersect  the  boundary  of  D. 
begin 

(1)  W  =  0;  1  =  0; 

repeat 

(2)  Compute  the  degree  2  parametric  approximant  $(s,i)  of 
F(x)  =  0  at  x; 

(3)  Compute  the  intersections  of  (s,t)  with  edges  of  B\ 

(4)  Refine  the  intersections  back  to  the  projected  surface  Sj  along  the 
edges  of  B; 

(5)  face{B)  =  {  [F,.<  (en,  pa ),  (etl,  (s,, .  i  >.<  I  el2.  p,2),  (e,2,  (s,2,  fl2 ) )  >]  | 

F,  is  a  face  of  B.  etl  a:u:  •  _  .ir^  ^uees  of  face  F,. 

Sf  intersects  ei;  at  p,  •••  p.  ;s  the  point  refined 
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from  and  is  on  etJ,  j  =  1,2  }; 

(6)  I  =  L  U  {<  B,  $(s,  t),  face(B)  >}; 

(7)  Remove  those  faces  in  face(B)  that  are  outside  £); 

(8)  for  B  €  W  begin 

(8-1)  if  B  and  B  have  a  common  face  F 

and  common  surface  intersections  px  and  p2 
then  begin 

(S'2)  face(B)  =  face{B)  -  {[/%<  (ex >, 

<  (e2,P2),(e2,(s2,t2))  >]}; 

(8-3)  face(B)  =  face(B)  -  {[F,<  (e1,px),{e1,{$1,tl))  >, 

<  (e2,ft),(cj,(5j,?2))  >]}; 

(8-4)  if  face(B)  =  0  then  W  =  W  -  {&}; 

end  /*  if  */ 
end  /*  for  */ 


(9) 

if  face(B)  7^  0 

then  begin 

(10) 

W  =  Wu{B}; 

(ID 

B  =  B: 

end 

(12) 

else  Select  a  cube  B  from  W: 

(13) 

Select  [F,  <  (ej 

,PiM«i.(«iifi))  >  <  (c2,P2),(c2,(32,t2))  >]  €  face(B); 

(14) 

Determine  ( u0 , 

u0)  from  (5i,ti)  and  (s2,  t2)  such  that  $(u0,  u o)  can  be 

refined  to  a  surface  point  x  =  (xt,x2)  where  X!  is  in  the  cube  B 
that  shares  F  with  B\ 
until  XV  =  0 


end 
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4.3  Detailed  Computations 

4.3.1  Intersecting  A  Parametric  Approximant  with  Edges  of  a  Cube 

A  crucial  step  in  the  Algorithm  4.1  is  the  computation  of  points  at  which  the 
parametric  approximant  intersects  the  cube’s  edges.  This  computation  is  more 
complicated  than  intersection  with  an  implicit  approximant.  With  an  implicit 
surface  approximant,  vertex  evaluation  can  be  used  to  determine  which  edges  in¬ 
tersect  the  surface  and  then  the  intersection  point  on  each  edge  can  be  estimated 
by  subdivision  or  interpolation  [17,  45,  Sj. 

When  tracing  the  intersections  of  $t(.s,/)  with  B ,  at  least  one  closed  loop  of 
curve  segments  will  be  formed  and  each  such  loop  corresponds  to  a  closed  loop  in 
the  (s,t)  parameter  space.  Since  $i(s,t)  approximates  the  surface  locally  within 
a  cube  B  and  the  intersection  points  along  the  edges  of  B  are  only  used  as  ap¬ 
proximates  to  surface  intersections  of  5/  with  edges  of  B ,  only  the  “nearest”  in¬ 
tersections  with  $i(s,f)  are  of  intertet.  More  precisely,  only  those  intersections 
whose  parameter  values  lie  on  the  innermost  loop  containing  (0,0)  are  considered. 
In  the  following,  when  we  mention  surface  intersections  with  edges  of  a  cube,  only 
surface-edge  intersections  of  this  kind  are  referred  to.  It  is  easy  to  see  that  $i(s.  t) 
may  intersect  the  edges  of  B  in  3,4,5,  or  G  points,  as  shown  in  Figure  4.1.  As 
degenerate  cases.  <I>i(s,t)  can  penetrate  a  face  of  B  without  intersecting  any  edge 
or  intersect  only  one  edge  of  B  in  two  points,  see  Figure  4.2.  We  formulate  the 
problem  of  finding  the  intersections  as  follows: 

Problem  4.1  For  given  $i(s,£)  =  (<i>i(s,  t),  <£3(3,  <))  °f  degree  2  and  a 

cube  B  containing  $i(Q.O),  compute  the  edge  intersection  points  of  $i(s.t)  with 
B. 

As  a  subproblem  of  (4.1).  we  explain  how  to  compute  the  intersection  of  $i(s,<) 
with  an  edge  of  B.  The  line  L  containing  an  edge  uu  of  the  cube  B  is  the  intersec¬ 
tion  of  two  face  planes  Pl(xi,x1,  x3)  =  0  and  P1(xl.xi,  x3)  =  0.  Thus,  substitution 
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(b):  2  points 


Figure  4.2  Intersecting  $x(s,t)  with  a  cube  B  (Degenerate  cases) 


of  $1(5,  t)  into  Pi  and  P2  yields  the  two  equations 


Pi(s,t)  =  0 
P2(M)  =  0 


(4.2) 


which  represent  the  intersections  of  ^(-M)  with  L  in  parameter  space.  For 
quadratic  approximants,  Pi(s,t)  and  p?(3,t)  are  of  degree  2.  System  (4.2)  can 
be  solved  either  algebraically  or  numerically. 

Using  a  resultant  method,  we  obtain  the  resultant  of  px  and  p2  as  a  degree  4 
polynomial  in,  e.g.,  s.  This  univariate  polynomial  is  then  solved  by  a  root  finding 
algorithm.  Substituting  each  root  s  into  px{s,t)  =  0  and  P2{s,t)  =  0  results  in  a 
polynomial  in  t  which  yields  solutions  t.  The  computed  solutions  (s,<)  are  checked 
to  see  if  $x(s,t)  is  on  uv.  Since  there  might  be  more  than  one  such  (3,  £),  it  is  not 
trivial  to  determine  which  (3,  t)  is  nearest.  The  computation  has  to  be  done  for  all 
edges  of  B  in  order  to  find  the  desired  intersections  of  $i(s,<)  along  edges. 

To  apply  Newton  iteration  to  system  (4.2).  an  initial  point  is  needed.  Possible 
initial  points  would  be  (s0,  to)  values  arising  as  the  intersection  of  the  tangent  plane 
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(a)  (b) 

a  :  intersections  of  tangent  plane  with  edges 
•  :  intersections  of  approximant  with  edges 


Figure  4.3  Intersecting  a  tangent  plane  with  the  cube  B 

of  ^(O,  0)  with  the  edges.  However,  these  initial  points  represented  as  (s0,f0)  >n 
the  parameter  space  might  not  be  close  enough  to  the  true  intersections  if  the  cube 
is  large  compared  to  the  radius  of  curvature  of  See  Figure  4.3(a).  Moreover, 
we  might  have  fewer  initial  points  than  the  number  of  intersections  with  $i.  and 
it  is  unclear  how  to  ascertain  that  all  intersections  have  been  found.  See  also 
Figure  4.3(b).  Notice  that 

1.  If  the  isoparametric  curves  $x(.s,0)  or  ^(O.t)  intersect  a  face  F  of  cube  B 
then  $i(s,f)  intersects  edges  of  F  at  least  at  2  points. 

2.  All  intersections  of  ^(s.f)  to  edges  of  B  form  a  closed  loop. 

With  the  above  observations,  we  derive  the  following  algorithm,  which  is  reliable 
and  fairly  efficient. 
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Algorithm  4.2  (For  Problem  4.1) 

Input 

B:  a  cube 

$x(s,t)  a  degree  2  parametric  surface  with  ^(O,!))  inside  B. 

begin 

(1)  C  =  0; 

(2)  Find  an  intersection  ,  ti )  with  an  edge  ex  of  5, 

where  e-i  =  F0C\  F\,  Fo  and  Fi  are  faces  of  B: 

(3)  LF  =  Fq\  Le  —  ej;  Ls  —  Sj;  Lt  —  ti ; 

(4)  repeat 

(5)  Find.the  second  intersection  <&t(s2,  t2)  on  an  edge  e2  of  Fi, 

where  e2  =  Fx  H  F2,  F2  is  a  face  of  B\ 

(6)  C  =  C  U  {<  ((F0,  Fi),ei,(s1,fl)),((F1,F2),e2,(s2,t2))  >}; 

(7)  Fo  =  Fi;  Fi  =  F2;  ei  =  e2;  si  =  s2;  tx  =  <2; 
until  F0  =  LF  and  (si,tt)  =  (Ls.Lt); 

end 

For  step  (2)  of  Algorithm  4.2,  we  do  the  following: 

Algorithm  4.3 

1.  Find  a  face  Fo  of  B  that  $i(a,0)  intersects,  say  at  (s, 0): 

(a)  Substitute  $i(s,0)  into  the  equation  of  plane  Fo  and  solve  the  resulting 
degree  2  univariate  polynomial  in  s. 

( b)  Take  the  real  root  i.  if  there  is  one,  that  is  closest  to  zero  and  determine 

is  inside  Fo.  If  so,  we  have  found  the  intersection. 

2.  (a)  Substitute  $i(s,  <)  into  the  equation  of  plane  F0  obtaining  a  plane  curve 

p(s<  t)  =0. 

(  b)  Trace  p{s.  t)  =  0  starting  at  (.s.0)  until  one  edge  e\  =  uT  of  F0  is  crossed 
over  at  (.s.  t). 


(c)  Apply  Newton  iteration  to  System  (4.2)  and  refine  (s.F)  to  a  solution 
($i,  fi)  of  (4.2).  $j(sj,  ti)  is  the  surface  point  on  ej. 

For  step  (4)  of  Algorithm  4.2,  a  computation  similar  to  step  (2)  of  Algorithm  4.2 
can  be  used.  Note  that  the  initial  tracing  direction  can  be  determined  easily, 
namely  the  one  that  is  pointing  to  the  inside  of  the  face  F\.  For  tracing  plane 
curves,  the  algorithm  and  implementation  of  [14]  is  used. 


4.3.2  Intersecting  A  Projected  Surface  with  Edges  of  a  Cube 

The  intersections  of  the  projected  surface  Sj  with  the  edges  of  a  cube  B  not 
only  determine  the  faces  that  are  transversal  to  Sf  but  also  provide  a  PLA  of  S/ 
within  the  cube.  Let  Pi(xj,X2,x3)  =  0  and  P2(x j,x2,x3)  =  0  be  the  two  face 
planes  whose  intersection  contains  an  edge  e  of  B.  A  surface  point  (xj.x2.x3)  on 
e,  if  there  is  one,  together  with  its  footpoint  (x4,  ...,x„)  on  the  basis  surface(s) 
satisfies  the  following  system  of  n  equations  in  n  variables 

/l(Xl,X2,...,Xn)  =  0 

fl(.X liX2,...,Xn)  =  0 


/„— 2(Xj,X2, - X  n)  =  0 

Pj(Xj.X2,X3)  =  0 

P2(Xj,X2,X3)  =  0 


(4.3) 


We  compute  the  surface  points  of  S /  on  edges  of  B  by  applying  Newton  iteration 
to  (4.3)  using  the  intersections  of  $j(s,t)  with  edges  of  B  as  initial  points.  That 
is,  when  $j(si,fj)  is  a  point  on  an  edge  ex  of  B.  <J>(si,fj)  e  Rn  serves  as  an 
approximation  of  a  zero  of  (4.3)  whose  natural  projection  to  (xj,  X2,  x3)-space.  say 
p,  is  the  refined  point  on  an  edge  e,  of  B.  In  short,  we  say  $j(sj,  tx )  on  ej  is  refined 
to  p  on  et.  An  approximation  <lv  'i-  ;i  can  be  refined  to  a  surface  point  on  the 
same  edge  or  to  a  surface  point  on  an  a.’.-.u  ent  I'dse.  Thus,  the  approximation  can 
be  refined  to  one  point  or  two  point '  -n-e.vn  in  Figures  4.4  and  4.5  respectively.  In 
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Figure  4.4  An  approximate  point  is  refined  to  a  surface  point 

the  latter  case,  an  additional  transversal  face  F  to  Sj  is  introduced.  Moreover,  two 
approximates  $i(si,tt)  and  $i(s2,f2)  might  be  refined  to  the  same  surface  point; 
see  Figure  4.6.  In  this  case,  the  face  F  is  not  transversal  to  5/.  We  describe  later 
how  to  predict  to  which  edge  an  approximant  is  refined.  Nevertheless  when  an 
approximant  is  refined  to  an  edge  e  using  (4.3),  P\  and  P2  are  the  plane  equations 
for  two  faces  whose  intersection  is  e. 

When  an  approximate  point  <0  is  on  an  edge  e\  of  B ,  it  is  assumed  that 

the  projected  surface  5/  will  intersect  t\  or  its  adjacent  edges  e2,e3,e4,e5,  and 
further  subdivision  is  necessary  if  this  is  not  the  case.  When  Newton  iteration  fails 

to  refine  4>i(si,  £i)  to  a  point  on  e,,  i  =  1 . 5.  we  subdivide  the  cube  B  into  eight 

equal-sized  cubes  in  order  to  acquire  better  approximations.  The  subdivision  will 
be  discussed  in  Section  4.3.6. 

Algorithm  4.4  accepts  as  input  a  cube  B.  its  list  of  k  transversal  faces  to  $i(s.  <)• 
and  the  associated  intersection  points  on  each  face.  The  list  of  transversal  faces  is 
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in  the  form 

{  [F 17  ®il  7  (^il7  1  )  5*7  ^i2l  ( -3*27  ^t’2  )  ''’j  I  *  I7  •  •  •  7  k} 

where  F,  is  a  transversal  face,  and  e,j  is  the  edge  of  B  on  which  $i(s,j,  lies,  for 
j  =  1,2.  The  algorithm  considers  each  transversal  face  F, ,  refines  ^1(3,!,  ttl )  and 
$i(s,2,  U 2)  and  determines  if  F,  is  transversal  to  S /.  When  the  two  approximates 
are  refined  to  two  different  surface  points  on  edges  of  F,,  the  face  F,  is  transversal 
to  5/.  If  the  two  approximates  are  refined  to  the  same  surface  point  then  F, 
is  not  transversal  to  5/  and  should  be  excluded.  Finally,  if  the  refined  point  of 
$i(s,i,  f.i)  is  different  from  that  of  $i(s(,_i)2,  <(,-1)2)  then  the  face  containing  the 
two  refined  points  is  transversal  to  5/  and  should  be  added  to  the  transversal- face 
list  of  Sf.  Note  that  $1(3,1,  f,i)  does  not  have  to  be  refined  if,  on  the  previous  face, 
-1)27  ^(»— 1)2 )  's  refined  to  a  point  lying  on  e(,_i)2.  If  this  is  not  the  case,  we 
see  if  $1(3,1,  <,i)  can  be  refined  to  the  edge  e  of  Ft  that  is  perpendicular  to  e,i  and 
e(,_i)2.  If  so,  the  face  F  containing  e  and  e(t_ i)2,  rather  than  F;,  is  transversal  to 
Sf]  otherwise  $i(s,!,  t,i)  must  be  refined  to  e(,_i)2;  see  Figure  4.7.  Moreover,  once 
$1(3,1,  Ui)  is  refined  to  a  surface  point  on  an  edge  of  a  face  F,  $1(3, 2,  t,2)  can  only 
be  refined  to  a  point  on  an  edge  of  F. 

We  describe  the  algorithm  in  pseudo  code  as  follows. 

Algorithm  4.4 

Input 

B:  a  cube 

{  [F,,  <  e,i,  (3,1,  t,i)  >,<  e,2,  (3,2,  t,2)  >]  |  i  =  i . A’}:  list  of  transversal  faces 

to  $1(3.  t)  and  intersections  on  edges. 

Output 

facet  B):  list  of  transversal  faces  to  5/  md  intorsectiors  in  the  form  of 
{  [F,,  <  ( e,i ,  p,i ),  ( e,i ,  (s,i ,  t,i ))  >.  <  (e,-2,  p,2  1.  .  t.,  > )  >  ]  !  i  =  l . /) 

begin 

(  I )  facet  B)  =  0: 
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Figure  4.7  The  refinement  of  ftl) 


(2)  i  =  0; 

repeat 

(3)  i  =  i  +  1; 

(4)  refine  <,i )  to  a  point  pn  on  an  edge  e,i  of  B\ 

(5)  refine  4>i(si2>  ^2)  to  a  point  pl2  °n  an  edge  e,2  of  B; 
until  pa  ^  Pi 2 

Let  F  be  the  face  containing  e,i  and  e,2 

(6)  face(B)  =  face(B)  U  {  [F,  r  (en,Pu),(eu,('S«i»<.i))  >. 

<  (e,2,Pi2),(e,2,(s,2,f,2))  >]  }; 

(7)  do  while  i  +  1  <  k  begin 

fS)  i  =  i+  1; 

(9)  if  e(,_i)2  =  e(l_l)2  or  (^(l_l)2  ^  s,i  and  /  t,i) 

then  begin  /*  is  refined  to  */ 

refine  $i(.s,2<  *,2 )  to  a  point  pl2  on  an  edge  e,2  of  B\ 
facet  B)  =  facet  B)U  {  [F„<  p(l_U2), (e,i, (stl,  <u ))  >. 

<  (e,2-P.2).(«.2,(5,2^.2))  >]  }: 


I  10) 

I  ID 
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(12)  end 
else  begin 

(13)  refine  $i($ii,  t,i)  to  a  point  pn  on  an  edge  e,x  of  B\ 

( 14)  if  e,i  =  Fi  n  Fe ,  where  Fe  D  F,_x  =  e(i_1)2, 

then  begin  /*  $x(aa,tu)  is  refined  to  2  different  points  */ 

(15)  face(B)  =  face(B)  U 

{  [Fe,<  [e{  t— 1)2>  P(t— 1)2)>  (®(t— 1)2>  (■*(» — 1)2»  ^(1—1)2))  >< 

<  ( eix, pn ),(e,i, (aa,<,i))  >j  }; 

(16)  refine  $i(sj2,tt2)  to  a  point  p,-2  on  an  edge  e;2  of  B\ 

(17)  face(B)  =  face(B)  U{  [F;,<  (e,i,pu),(e«i,(s. >, 

<  (e,2,P,2),(e,'2,(s,2,t,2))  >]  }; 

end 

else  begin 

(18)  refine  $i(.Si2,tl2)  to  a  point  pt2  on  an  edge  et2  of  B\ 

(19)  if  pa  ^  pi2  then  begin 

(20)  let  F  be  the  face  containing  e,x  and  el2; 

(20)  face(B)  =  face(B)  u{  [F,<  (ea,pa),(e,i,(s,x,t,i))  >. 

<'  (®»2>  P«2  )<  (^i2»  (fl«2»  ^»2))  >  ]  }> 

end  /*then*/ 
end  /*else*/ 
end  /*else*/ 
end  /*do*/ 

end 

4.3.3  Strategies  for  Stepping 

After  transversal  faces  of  a  cube  B  have  been  found,  we  track  the  surface  5/  by 
propagating  cubes  along  transversal  edges,  taking  care  that  no  cube  is  processed 
twice.  When  considering  an  adjacent  transversal  cube  B.  a  surface  point  (xi.x2) 
such  that  Xj  6  B  is  required  for  computing  the  local  approximant  $(s.t)  from 
which  <J>i(s,<)  serves  as  the  approximation  of  Sj  inside  B.  To  find  such  a  surface 
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point  (xi,x2),  one  possibility  is  to  extend  the  surface  approximant  $i(s,t)  of  5/ 
in  B  to  B  and  locate  an  appropriate  (uo,  Vo)  such  that  $(u0,v 0)  can  be  refined  to 
the  desired  surface  point.  Since  the  local  shape  of  5/  within  a  cube  B  is  approx¬ 
imated  by  it  is  clear  that  the  location  of  Xi  closely  influences  the  overall 

performance  of  the  approximation.  Ideally,  we  expect  that  xt  is  approximately 
at  the  center  of  all  surface  intersections  on  edges  of  B.  The  approach  we  take  is 
to  choose  (u0,  v0)  so  that  {u0,v0)  is  the  barycenter  of  (ux,  vx ),...,  (u*,  u*,)  where 
$i(u j,  vx ), . . . ,  l>i (ufc,  Vk)  are  intersections  of  $i(s,  t )  with  edges  of  B. 

Recall  that  with  each  transversal  face  F  of  B  there  is  an  associated  record 

[F,  <  (et,pi),(e1,(sll  ti))  >.  <  (e2,f  [e2,  (s2 ,h))  >] 

where,  for  i  =  1,2,  p,  is  the  intersection  of  5/  on  edge  e,  of  F,  p,  is  refined  from 
$i(s,-,  ti),  and  e,  is  the  edge  of  B  with  which  $i(s,  t)  intersects  at  $i(s,,  t,).  When 
e!  or  e2  is  not  on  F,  it  is  clear  that  is  no  longer  appropriate  for  the  above 
computation.  In  this  case,  we  first  compute  s0  =  j(si  +  s2)  and  t0  —  %(tx  +  t2) 
and  refine  ^(so,  ^o)  to  a  surface  point  p  of  5/  on  face  F  and  then  on  p  we  derive 
a  new  approximant  $i(s,£);  see  Figure  4.8(a).  This  computation  is  also  applied 
to  the  case  in  which  (si,  £j)  =  (s2,  t2);  see  also  Figure  4.8(b).  Moreover,  for  better 
approximation  of  ^(s.t)  within  B.  this  computation  can  generally  be  applied 
when  both  e\  and  e2  are  edges  of  the  face  F.  Hence,  as  an  uniform  approach  for 
all  cases,  we  compute  such  a  new  approximant  <&x(s,t)  for  the  transversal  face  F. 

Nevertheless,  Newton  iteration  may  fail  to  refine  the  computed  $(u0,  v0)  to  a 
surface  point  (x1(x2)  or  the  iteration  succeeds  but  Xi  is  outside  B.  When  this 
happens,  the  following  heuristic  approach  can  be  applied: 

1.  compute  s0  =  j(so  +  u0)  and  t0  =  j( tQ  +  v0). 

2.  refine  4>i(s0,  t0)  to  a  surface  point  on  which  a  new  approximant  $i(s.£)  is 
derived. 


3.  compute  (  u0,  r0). 
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4.  refine  $i(uoi  vo)  to  a  surface  point. 

5.  go  to  step  1  if  the  step  4  fails. 

4.3.4  Newton  Iteration 

Newton  iteration  has  been  used  in  this  chapter  to  refine  an  approximate  zero 
p0  to  a  true  zero  p  of  the  system  of  equations  by  generating  a  sequence  of  points 
Pi,  p2, . . . ,  — p.  Systems  (4.2)  and  (4.3)  are  a  0-dimensional  nonlinear  systems 
G(x)  =  0,  where  G  :  Rm  — +  Rm  for  some  m.  The  Newton  iteration  for  solving 
such  systems  is  given  by 

Z?G(pfc)(pi+1  -  pfc)  =  -G(pfc) 

This  is  a  linear  system  of  equations  for  Pfc+i,  and  if  DG(pk)  is  nonsingular,  it  can 
be  solved  by  general  linear  solvers. 

When  refining  an  approximate  to  a  surface  point  of  F(x)  =  0  by  Newton 
iteration,  we  solve  the  following  underdetermined  linear  system 

V/,(pfc)  •  Afc  = -/,(pfc),  *  =  1.2 . n  —  2  (4.4) 

where  A/t  =  pt+i  —  p/t  and  DF(pfc)  is  of  dimension  (n  —  2)  x  n.  Equation  (4.4)  is 

the  same  as  equation  (3.12).  When  D F(pfc)  has  rank  n  —  2.  the  general  solution  is 

Afc  =  Oriti  +  0^2  ~f  +  •  •  •  +  Pfc )  (4.5) 

where  ti  and  t2  form  a  unit  vector  base  of  the  tangent  space  of  F(x)  =0  at  p*. 

If  it  is  wished  to  make  the  computation  more  numerically  stable,  one  may 
use  singular  value  decomposition  as  we  have  done  in  Section  3.4.  The  matrix 
I'DFfp*))7'  is  factored  as  ( Z?F(pfc ) ) r  =  CEVT,  where  U  =  [U|,...,un]  €  R',xn 
and  V  =  [v1? . . . ,  vn_2]  €  R(n~2>x<n-2)  are  orthogonal  matrices,  and  E  €  Rn*in-2) 
is  a  diagonal  matrix.  Directly  substituting  the  factorization  to  system  (4.4)  we 
obtain 


VZrrr\,  =  -Ffp;.) 


whose  solution  can  be  generally  written  as 


Afc  =  7iUi  +  72U2  - h  7„u„ 


where  7i,...,7„_2  are  uniquely  determined  by  7,  =  (-vf F(pi))/E^  and  7,,., 
and  7n  are  arbitrary.  For  A*,  we  assign  7n-i  =  7„  =  0  so  that  A*  is  in  the 
normal  space  of  Sp  at  p*  since  Ui, . . . ,  un_2  span  the  same  space  as  the  gradients 
WiPfc)i  •  •  •  1  y/n-2(Pfc)-  Alternatively,  we  can  compute  the  QR-factorization  of 
(DF(pfc))T  as 

(DF(pfc))T  =  Q 

where  Q  is  a  n  x  n  orthogonal  matrix  and  R  is  a  (n  —  2)  x  (n  —  2)  nonsingular, 
upper  triangular  matrix.  To  compute  Afc,  we 

1.  solve  R7y  =  -F(pfc)  for  y  6  Rn*2. 

2.  solve  QrAfc  =  (y,0)r,  i.e.,  A*  =  Q(y,0)r. 

Notice  that  the  last  2  columns  of  Q  forms  an  orthogonal  basis  of  the  tangent  space 
of  5f*  at  pfc.  Thus  A*  so  computed  is  a  linear  combination  of  the  gradients. 

4.3.5  Some  Error  Analyses 

A  problem  of  practical  and  theoretical  interest  is  to  estimate  the  error  bound 
between  the  projected  surface  and  its  local  parametric  approximant  within  a  cube 
or  on  cube’s  edges.  The  computation  of  this  estimate  would  depend  on  error 
analyses  of  both  the  projection  and  the  approximation.  It  is  not  clear  now  how 
to  approach  the  above  problem.  Instead,  we  consider  the  estimation  of  the  error 
bound  between  the  local  parameterization  of  the  projected  surface  and  the  local 
parametric  approximant  within  a  cube,  under  the  assumption  that  the  neighbor¬ 
hood  of  convergence  of  the  local  parameterization  covers  the  parametric  points 


R 
.  0 
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Recall  that  the  local  parametric  approximant  $1(3,  t)  of  Sf  consists  of  the  initial 
terms  of  the  Taylor  series  of  the  local  parameterization  up  to  degree  2.  Let  *i(.s,  t) 
denote  the  local  parameterization  of  Sf.  Consider  the  coordinate  functions  0,  (•M) 
and  ipi(s,t)  of  *1(3,  £)  and  *i(s,f)  respectively.  From  the  differential  calculus. 

0i(s,  t)  =  ©i(s,  t)  +  R\3)(rt,  s,  t) 


for  some  0  <  r,  <  1,  where  R[J\ri,s,  t)  =  jr((s,t)  ■  (Jj,  £)yii\(rt(s,t)).  Thus 
IvdM)  -  Oi(s.  f)|  =  j R[3)(T,,s,t)\  <  max  \R\3){t„  s.  <)| 

0<r,<l 

However,  0,(s,t)  is  unknown  and  hence  the  above  bound  cannot  be  computed. 
Consider  one  exact  representation  of  the  error 


v,{s,t)  -  Ot(s,t)  =  T,(3'(s,t)  +  R(t4)(Tt,sJ) 

for  0  <  r,  <  1,  where  TjJ\s,t)  =  jr((s,t)  •  ( iiW'i’iis,  0-  With  h  =  ||(s,f)||, 
T-3\s,  t)  is  an  0{h 4)  accurate  approximation  of  the  error  v>,(s,  t)~  <t>,(s,  f)-  Thus  the 
vector  (T{3)(s.  t).  r23)(s,  t ),  Tj3)(s,  t))T  estimates  *i(s,  t)  — *1(3,  t)y  componentwise, 
to  within  0(h*)  accuracy,  and  hence 

||*1(s,t)-$1(s,t)||  =  i|(T1(3,(s,o.r‘3)(5.t).T‘3)(5,t))r||  +  0(^) 

/  o  \ 

Notice  that  the  coefficients  of  T,  (s.t)  can  be  computed  using  equation  (3.12). 


4.3.6  Adaptive  Subdivision 

The  cube  size  6  is  an  input  to  the  PLA  algorithm.  To  determine  an  appropriate 

5  for  a  given  problem  a-pnorx  is  nontrivial.  In  practice,  we  would  expect  a  large 

6  initially,  and  perform  adaptive  subdivision  whenever  necessary.  There  are  cases 
in  which  a  cube  must  be  subdivided  in  order  to  continue  the  computation.  The 
crises  are  as  follows. 

1.  When  a  parametric  approxim.mt  ;>»*ner rates  a  face  without  intersecting  any 
edge  of  the  face  as  shown  in  F:r:re  12. 
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q  qi 


□:  edge  intersections  of  the  approximant. 
•  :  edge  intersections  of  5/. 

Pi  and  p2  are  refined  to  p 
7i  and  q?  are  refined  to  q 


An  invalid  polygon 


2.  When  Newton  iteration  fails  to  refine  the  intersections  of  the  parametric 
approximant  with  edges  of  the  cube  to  surface  points  on  edges. 

3.  When  an  approximate  point  is  refined  to  an  unexpected  surface  point  and 
hence  an  invalid  polygon  is  produced,  e.g.,  see  Figure  4.9. 

Like  PL  A  methods  that  use  vertex  evaluation,  some  portions  of  the  surface  might 
be  truncated  if  we  do  not  subdivide.  Two  typical  examples  are 

1.  The  projected  surface  5/  penetrates  a  face  in  a  closed  curve  without  inter¬ 
secting  the  boundary  edges;  see  Figure  4.2(a).  In  this  case,  both  neighboring 
cubes  sharing  that  face  should  be  subdivided.  This  is  difficult  to  detect  in 
general. 

2.  5/  intersects  only  one  edge  of  a  face:  i.e.,  Sj  intersects  an  edge  of  a  face 
ut  two  points  and  has  no  intersections  with  others:  e.g.,  see  Figure  4.2(b). 
Thus  all  four  neighboring  cubes  along  that  edge  are  subdivided  for  a  reliable 
polvgonization. 
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Figure  4.10  A  crack  on  the  shared  face 

Subdivision  is  performed  by  calling  Algorithm  4.1  recursively,  with  the  parent 
cube  as  the  domain  of  interest  and  with  the  subcube  containing  the  regular  surface 
point  as  the  starting  cube.  Because  of  size  differences,  cracks  may  occur  along  the 
transversal  faces  of  the  parent  cube.  See  also  Figure  4.10.  In  order  to  close  such 
cracks,  a  data  structure  is  needed  that  records  topological  adjacencies  of  cubes. 

We  modify  the  data  structures  of  Algorithm  4.1  by  associating  with  each  parent 
cube  a  List  of  transversal  subcubes,  and  pointers  to  neighboring  transversal  cubes 
at  the  same  level.  Thus  for  each  cube  B  we  have 

<  B,$(s,t),face(B),(pti,ptt,pt3,pt4,pt$,pt6),ptp,LB  >  (4.6) 

where  ptv  is  a  pointer  to  the  parent  cube,  Lg,  a.  list  of  elements  in  the  form  (4.6). 
is  the  list  of  transversal  subcubes  of  B ,  and  pti  is  a  pointer  to  the  adjacent  cube 
sharing  face  F,  with  B.  The  list  Lg  is  assigned  to  the  parent  cube  on  the  return 
from  subdivision.  The  pointers  pt,  to  neighboring  transversal  cubes  are  assigned 
when  such  adjacencies  are  confirmed  in  step  S  of  Algorithm  4.1.  The  face(B)  is 
derived  from  Lg  if  B  is  subdivided. 
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4.3.7  Polygonization  and  Local  Refinement 

The  PL  A  of  5/  is  a  polygonal  representation  given  as  a  list  of  polygons  derived 
from  those  cubes  that  intersect  Sf.  A  list  of  transversal  cubes  lying  in  a  given 
domain  D  is  produced  by  Algorithm  4.1.  For  each  transversal  cube,  say  B.  there 
is  an  associated  list  of  transversal  faces  and  surface  intersections  from  which  the 
polygon  approximating  5/  within  B  is  formed. 

When  adaptive  subdivision  is  incorporated  into  Algorithm  4.1.  a  transversal 
cube  B  might  be  recursively  subdivided  and  hence  the  polygonization  processing 
is  confined  to  terminal  transversal  subcubes.  As  mentioned  in  Section  4.3.6.  the 
polygonization  scheme  must  close  cracks  along  the  shared  face  of  two  transversal 
cubes. 

Recall  that  we  associated  with  a  face  F  of  a  cube  C  a  record  of  intersection 
information  in  the  following  form 

[F,<  (ej,pj),  (ej,  (sj,  ti))  >,<  (e2, P2),  (e2,  (s2,  t2))  >} 

When  tracing  the  transversal  faces  of  a  cube  C  in  order  to  form  a  polygon  over 
C.  we  check  on  each  transversal  face  F  to  see  if  the  adjacent  cube  C  has  been 
subdivided.  If  C  is  subdivided,  possibly  more  than  once,  we  trace  the  surface-edge 
intersections  on  face  F  of  C  starting  from  px  on  ex.  If  the  faces  of  a  cube  are 
numbered  consistently,  this  computation  can  be  achieved  efficiently  by  processing 
the  faces  of  the  same  number  on  the  subcubes  of  Lg  that  share  F  with  C .  Note 
that  when  C  is  a  subcube  of  B  while  C  is  not.  C  can  be  located  by  following  the 
parent  pointers  and  the  pointers  to  the  neighboring  cubes. 

Once  a  polygon  is  formed,  it  is  desirable  to  refine  locally  the  polygon  according 
to  some  criterion  provided  by  the  user.  e.g..  the  maximum  deviation  of  vertex 
normals  from  the  normal  of  p  =  <^(0,0).  L>*r  the  polygon  P  be  the  list  of  vertices 

(pi,p2 . Pfc j,  for  some  k.  and  let  A,  be  the  imt  normal  of  Sf  at  p,,  1  =  1 . k. 

Also  let  .V  be  the  unit  normal  at  p.  As  <h>>wn  .:i 


he  maximum  deviation  of 
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vertex  normals  from  the  normal  of  the  p  can  be  estimated  by 


max(iV,  •  N) 
i<'<* 

When  the  maximum  deviation  exceeds  some  tolerance,  polygon  P  is  replaced  by 
triangles  [p,pi,pi+ ij  for  i  =  1 _ ,  k  —  1.  and  each  of  the  triangles  is  refined  ac¬ 

cordingly.  Note  that  on  the  refinement  of  triangles,  the  midpoint  of  the  triangle  is 
refined  to  a  surface  point  by  Newton  iteration. 

Surface  normal  computations  are  also  needed  when  shading  the  surface  to  pro¬ 
duce  an  image.  For  computing  the  normals  of  the  projected  surface  5/  directly 
from  F(x)  =  0,  see  Section  3.1.1.  Notice  that  the  coordinates  of  the  polygon  ver¬ 
tices  in  Rn  are  stored  in  order  to  do  the  refinement  and  compute  surface  normals. 
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5.  CONCLUSIONS  AND  FUTURE  WORK 

5.1  Local  Implicit  Approximations  of  Curves  and  Surfaces 

We  have  presented  a  method  for  computing  a  local  implicit  approximation 
of  parametric  curves  and  surfaces.  The  method  works  for  both  polvnomially  and 
rationally  parameterized  curves  and  surfaces,  and  achieves  an  order  of  contact  that 
can  be  prescribed.  In  the  case  of  nonsingular  curve  points,  the  approximant  must 
be  irreducible,  but  in  the  surface  case  additional  safeguards  have  been  incorporated 
into  the  algorithm  to  ensure  irreducibility.  The  method  also  yields  meaningful 
results  for  many  types  of  singularity.  The  algorithm  is  capable  of  determining  the 
exact  implicit  form  without  extraneous  factors  when  the  approximant  is  formulated 
with  the  exact  degree  of  the  implicit  form. 

The  method  provides  a  middle  ground  between  two  major  approaches  for  eval¬ 
uating  the  intersection  curve  of  two  parametric  surfaces,  that  is,  subdivision  and 
substitution  methods.  It  is  well  known  that  subdivision  methods  are  robust  and 
can  locate  all  intersection  branches,  but  at  the  expense  of  creating  a  large  vol¬ 
ume  of  data,  while  the  substitution  method  provides  an  exact  representation  of 
the  intersection  with  the  help  of  often  expensive  implicitization  techniques.  In 
the  context  of  subdivision  methods,  our  implicit  approximations  have  the  poten¬ 
tial  of  reducing  the  number  of  generated  surface  approximants  since  we  are  not 
restricted  to  only  planar  approximants.  In  the  context  of  substitution  methods, 
the  approximations  avoid  the  high  cost  of  implicitization.  In  both  cases  a  number 
of  practical  issues  remain  open  for  exploration,  including  the  trade-off  between 
the  degree  of  the  approximant  and  the  accuracy  with  which  the  curve  or  surface 
has  been  approximated.  In  particular,  a  comparative  evaluation  of  our  method  is 
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desirable  that  contrasts  its  performance  with  other  surface  intersection  methods, 
such  as  the  one  based  on  subdivision. 

5.2  Local  Approximations  of  2-Surfaces 

We  have  demonstrated  that  certain  surfaces,  including  offsets,  blends  and 
Voronoi  surfaces,  can  be  formulated  as  the  projection,  into  3-space,  of  2-surfaces 
F(x)  =  0  in  Rn.  For  such  surfaces,  we  have  proposed  several  computation  schemes 
that  (1)  describe  the  local  geometry,  including  computations  of  normal  vectors, 
tangent  vectors,  and  normal  curvatures,  of  the  projected  surface,  and  (2)  derive 
degree  two  local  implicit,  local  explicit  and  local  parametric  approximations  of 
the  projected  surface.  We  believe  that  these  methods  will  be  of  practical  interest 
in  rendering  2-surfaces  in  high-dimensional  space  and  computing  surface/surface 
intersection. 

In  computing  the  degree  two  implicit  approximation  of  the  projected  surface,  an 
3x10  linear  system  is  obtained,  which  leaves  us  two  degrees  of  freedom.  Problems 
of  practical  and  theoretical  interest  include  describing  the  set  of  approximants 
parameterized  by  the  two  degrees  of  freedom  and  deriving  criteria  from  which  to 
determine  a  member  of  the  family  of  approximants  to  select. 

Several  other  problems  are  of  interest,  for  example,  computing  the  Gaussian 
and  mean  curvatures  of  the  projected  surface  at  a  point  directly  from  F(x)  =  0 
without  variable  elimination,  and  computing  parametric  approximations  of  the 
projected  surface  at  singular  surface  points. 

5.3  Piecewise  Approximations  of  2-Surfaces 

The  ability  to  derive  a  piecewise  linear  approximation  (PLA)  of  a  2-surface 
defined  implicitly  in  Rn  should  be  essential  in  the  interactive  design  environment 
since  the  PLA  allows  one  to  take  advantage  of  hardware  capabilities  and  reduces 
the  cost  of  expensive  ray  tracing  in  the  rendering  process. 
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We  have  presented  an  algorithm  that  computes  the  PLA  of  a  2-surface  defined 
bv  a  system  of  nonlinear  equations  in  n  variables,  where  n  >  3,  but  whose  natural 
projection  to  3-space  is  the  surface  of  interest.  This  is  an  algorithm  that  deals 
with  a  2-surface  in  Rn,  but  performs  all  major  computations  in  3-space.  Hence  its 
performance  on  computing  the  PLA  of  complex  surfaces,  including  offsets,  blends, 
and  Voronoi  surfaces,  is  much  more  efficient  than  methods  that  work  in  n-space. 

A  number  of  issues  await  future  research.  Our  algorithm  takes  as  input  a  cube 
of  a  prescribed  size  containing  a  given  point  on  the  (projected)  surface.  One  might 
ask  how  one  should  determine  the  cube  size.  In  the  process  of  the  algorithm,  a 
degree  two  parametric  approximant  is  derived  to  approximate  the  projected  surface 
within  the  cube  and  its  intersections  with  the  cube's  edges  serve  as  initial  points 
to  be  refined  to  the  intersections  of  the  projected  surface  with  edges  of  the  cube. 
Moreover,  cube  subdivision  can  be  easily  applied  when  necessary  and  due  to  the 
availability  of  the  parametric  approximant  the  local  refinement  within  a  polygon 
can  be  efficiently  performed.  Thus,  we  would  expect  a  large  cube  size  initially, 
followed  by  adaptive  subdivision.  An  ideal  cube  size  would  balance  local  geometry 
and  global  geometry  such  that  only  a  small  number  of  cube  subdivisions  are  needed. 
Such  trade-offs  remain  an  important  problem  to  be  explored  in  greater  depth. 

Locating  a  seed  point  on  the  projected  surface  is  in  general  difficult  especially 
for  blends  and  Voronoi  surfaces.  Space  decomposition  is  useful  for  surfaces  in  R3, 
but  seems  expensive  for  2-surfaces  in  R"  where  n  is  large.  It  remains  an  urgent 
challenge  for  PLA  algorithm  based  on  continuation  methods. 

An  efficient  and  reliable  way  to  estimate  stepping  into  an  adjacent  transversal 
cube  is  crucial  to  the  performance  of  our  algorithm.  Although  we  have  given 
a  heuristic  approach  that  is  based  on  the  local  parametric  approximant.  a  more 
quantitative  method  would  be  desirable  that  accounts  for  the  local  geometry  of 
the  projected  surface. 
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